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ABSTRACT 


The  analysis  of  a fly-by-wire  longitudinal  control 
system,  specifically  that  of  the  space  shuttle  or'oiter, 
was  undertaken  in  order  to  demonstrate  the  construction 
of  a mathematical  model  depicting  the  relationships 

forcing  ‘function.  3.r.d  • 5l3.cn  "fsccc  of* 

modern  control  theory,  including  stability,  was  devel- 
oped. Several  computer  programs  were  written,  which 
should  be  of  value  to  the  Department  of  Aeronautics, 
for  the  HP9330  comou ter/plot ter ; these  programs  are  basic 
to  the  study  of  control  theory,  demonstrate  the  importance 
of  the  transfer  function,  the  characteristic  equation,  and 
the  various  forms  of  feedback,  and  will  plot  time  and 
frequency  (Bode)  response ' graphs  saver,  the  proper  inputs. 
The  Continuous  System  Modeling  Program,  version  III,  ar.d 
the  IBr’360  were  used  to  analyze  the  complex  control  system 
installed  in  the  Crbiter.  The  demonstration  of  the  model 
and  its  interface  with  the  C3TP  program,  was  given,  and  the 
efficiency  of  this  procedure  was  made  clear. 
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I.  PHYSICAL  DESCRIPTION  OF  THE  OR8ITSR 


The  space  shuttle  orbiter  is  the  vanguard  of  a new 
generation  and  new  concept  of  spacecraft,  having  been 
designed  to  fly  as  an  atmospheric  vehicle  as  well  as  to 
operate  as  a spacecraft  in  near-earth  orbit.  In  order  to 
function  in  this  dual  role,  it  requires  the  performance, 
stability,  and  control  systems  of  both  types. 

For  orbital  flight,  attitude  control  is  achieved  by 
means  of  a reaction  control  system  such  as  those  used  in 
earlier  spacecraft.  Whereas  earlier  vehicles  followed  an 
unguided  re-entry  and  parachute  let-down,  the  orbiter  is  to 
be  flown  to  a runway  landing  and  therefore  requires  the 
conventional  aerodynamic  roll,  pitch,  and  yaw  controls 
for  atmospheric  flight.  The  longitudinal  system  for 
maintaining  the  desired  pitch  angle,  pitch  rate,  and 
stability  about  the  pitch  axis  is  the  subject  of  this 
paper. 

The  space  shuttle  orbiter  is  of  conventional  design, 
incorporating-  a low-mounted,  highly  swept  delta  wing  and 
a small  vertical  stabilizer  on  a thick,  rather  square  body. 
Sea  level  weight  of  the  vehicle,  in  re-entry  configuration, 
is  approximately  186,000  lbs.  The  wetted  area  of  the  wing 
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is  2690  square  feet,  yielding  a wing  loading-  of  69  psf. 
Overall  length  is  122.3  ft,  and  wing  span  is  78  ft.1 
The  landing  speed  of  the  orbiter,  necessarily  high  because 
of  the'  wing  loading,  is  190  knots. 

The  primary  aerodynamic  control  surfaces  are  split 
elevons  at  the  trailing  edge  (to  perform  the  functions  of 
ailerons  and  elevators)  and  a combination  rudder/speed- 
brake  at  the  trailing  edge  of  the  vertical  stabilizer. 

The  rudder  splits  symmetrically  along  its  vertical  hinge 
to  produce  aerodynamic  drag  for  speed  control,  hence  its 
use  as  a speed  brake.  Body  flaps  are  located  at  the 
trailing  edge  of  the  fuselage  and,  although  they  were 
originally  intended  to  provide  heat  shielding  for  the 
main  engines  during  re-entry,  will  be  used  as  an  active 
longitudinal  trim  device.  Figure  1-1  is  a three-view  of 
the  vehicle  and  shows  the  location  of  the  primary  control 
surfaces. 

The  elevons  are  split  because  of  wing  deflection  and 
hinge  moment  considerations,  and  are  quite  large,  comprising 
twenty  per  cent  of  the  total  wing  area.  Their  maximum  de- 
flection is  35  degrees  up  and  20  down;  maximum  deflection 

1 By  comparison,  the  commercial  Douglas  DC-9  measures 
119.3  feet  in  length  and  has  a span  of  93»^  feet. 
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Figure  1-1 
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rate  is  twenty  degrees  per  second. 

The  speed  brake  provides  additional  longitudinal 
control  in  some  regimes  by  generating  a pitching  moment, 
but  its  primary  use  is  for  terminal  area  energy  management. 

The  elevons  are  the  primary  pitch  attitude  control, 
but  since  their  normal  deflection  is  trailing  edge  down, 
the  body  flap  is  available  to  reduce  the  load  on  them. 

Use  of  the  flap  enables  the  elevons  to  trail  as  much  as 
possible  and  minimizes  their  chances  of  overheating. 

The  body  flap  is  also  an  aid  to  longitudinal  control  at 
high  angles  of  attack  where  the  elevons  lose  much  of 
their  effectiveness.  Due  to  rudder  blanking,  the  orbiter 
is  in  fact  directionally  unstable  in  this  regime  (the 
early  re-entry  phase).  The  existence  of  a lightly  damped 
dutch  roll  mode,  however,  provides  dynamic  stability  for 
the  vehicle. 

The  trajectory  followed  by  the  vehicle  from  orbit 
to  touchdown  is  sensitive  and  complex  as  with  former 
spacecraft  and  occurs  as  follows.  Following  de-orbit, 
the  reaction  control  system  engines  orient  the  orbiter 
to  a normal  front-end-f orward  top-side-up  attitude  prior 
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to  reaching  re-entry  interface  at  400,000  feet.  The  angle 
of  attack  during  re-entry  varies  between  thirty  and  forty 
degrees.  When  the  velocity  decreases  to  900 0 feet  per 
second,  a slow  pitch  down  maneuver  commences  which  is  com- 
pleted by  70,000  feet;  here,  the  angle  of  attack  has  been 
reduced  to  ten  degrees  and  the  velocity  to  1500  feet  per 
second.  Terminal  area  energy  management  is  then  enabled 
to  control  airspeed  and  to  position  the  orbiter  to  arrive 
on  the  final  approach  at  12,000  feet,  six  nautical  miles 
from  touchdown.  Most  of  the  phase  prior  to  the  final 
approach  position  will  be  flown  at  a large  roll  angle  and 
at  relatively  high  "g"  in  order  to  follow  a steeper  tra- 
jectory. This  results  in  shortening  the  range  and  time  to 
touchdown,  which  reduces  the  total  heat  load  endured  by  the 
vehicle  and  extends  the  lifetime  of  the  thermal  croteccicr. 
system. 

The  flight  control  system  itself  has  three  modes: 
automatic,  control  stick  steering,  and  direct.  Each  mode 
may  be  selected  individually  for  longitudinal  and/or 

lateral-directional  control. 

The  automatic  mode  consists  of  hands-off  coupling 
of  guidance  commands  through  the  digital  flight  control 
system  to  the  control  surfaces.  Control  stick  steering 


is  a pitch  (roll)  rate  command  system  for  the  longitudinal 
(lateral-directional)  channel  with  gain  scheduled  as  a 
function  of  the  dynamic  pressure.  This  mode  is  further 
compensated  for  angle  of  bank  and  for  the  pitching:  moment 
created  by  the  speed  brake.  The  direct  system  relates 
controller  deflections  proportionally  to  the  control 
surfaces;  this  mode  will  be  hopefully  unnecessary  as  the 
vehicle  is  marginally  stable  longitudinally,  even  with  a 
forward  center  of  gravity.  "The  aerodynamic  characteristics 
the  orbiter,  together  with  its  total  dependence  on 
avionics  for  flight  control,  classify  it  as  a digital 
fly-by-wire,  control-configured  vehicle. The  entire 
approach,  from  de-orbit  to  re-entry  to  touchdown,  will 
normally  be  flown  entirely  in  the  automatic  mode. 


^ LCOL  C.O.  Fullerton,  USAF,  "Space  Shuttle  Orbiter  Approach 
and  Landing-  Test  Program,"  The  Society  of  Experimental 
Test  Pilots,  1 Report . page  1 53 
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LATERAL  ACCEL  Figure  1-2.  Block  diagram 

of  the  control  stick 
steering  mode 


II.  DEFINITIONS  AND  SYMBOLOGY 


As  in  any  field  of  study,  that  of  automatic  control 
has  many  terms  which  have  a particular  meaning  when  used 
in  this  context.  So  that  further  references  to  some  of 
these  terms  is  facilitated,  the  following-  definitions 
apply : 

Angle  of  attack:  the  angle  between  the  fuselage 
reference  line  and  the  flight  path;  denoted  by  oc . 

Block  diagram:  a graphical  representation,  such  as 
a schematic  diagram,  of  the  physical  components  of  a system 
or  a set  of  mathematical  equations  describing  their  input/ 
output  relationships. 

Bode  plot:  a plot  of  the  magnitude  of  response  of  a 
system  to  a sinusoidal  input  measured  against  the  fre- 
quency of  that  input. 

Characteristic  equation:  the  denominator  of  a system's 
transfer  function  set  equal  to  zero. 

Closed-loop:  a system  in  which  the  input  is  affected 
or  modified  by  some  stage  of  output. 

Control  matrix:  describes  how  the  control  incut  affects 
the  plant. 

Control  system:  a set  of  physical  components  related 


in  such  a manner  so  as  to  control  or  regulate  itself  or 


some  other  system. 

Damping:  frictional  retardation  of  the  magnitude  of 
a deflection  (displacement)  or  rate. 

Damping-  ratios  the  ratio  of  the  amount  of  damping 
actually  present  in  a system  to  the  amount  that  would 
critically  damp  that  system. 

Frror  signal:  the  actual  actuating  signal  of  the  system 
it  is  generally  the  difference  between  the  input  and  feed- 
back signals. 

Feedback:  a signal,  which  is  some  form  of  the  output 
signal,  that  is  added  to  or  subtracted  from  the  original 
input  in  order  to  compare  the  desired  output  with  that 
actually  obtained. 

Forcing  function:  see  Input. 

Frequency  response:  a measure  of  the  output  of  the 
system  as  a function  of  input  frequency. 

Dain:  the  measure  of  magnitude  of  an  input  or  feed- 
back signal. 

Initial  condition:  one  of  the  two  types  of  inputs  to 
a system. 

Input:  a stimulus  applied  to  a control  system  in 
order  to  produce  a desired  output. 

Linear  differential  equation:  one  consisting  of  the 
sum  of  linear  (first  degree)  terms;  it  has  the  property 
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that  the  response  of  a system  due  to  several  inputs  acting 
at  once  is  equal  to  the  sum  of  the  responses  of  each  input 
acting  alone.  An  equation  may  be  linearized  by  making 
valid  assumptions  (approximations)  about  the  system  or 
by  considering-  the  system  to  operate  in  a restricted 
environment  where  the  variations  from  linearity  are 
small  (piece-wise  linearity). 


Oven 

-loon:  a sys 

tern 

in  which  the 

input  is  not 

af  f ec  ted 

or  modified 

by 

any  stage  of 

ou  tout . 

Outout:  the  actual  response  of  the  system  caused 
by  the  input  reacting  with  the  various  components  of  the 
system;  it  is  not  necessarily  the  response  implied  by  the 
input. 

Flant:  the  object  or  system  which  is  to  be  controlled. 

Reaction  control  system:  a control  system  which 
generates  rolling,  pitching,  and/or  yawing  moments  by 
venting  jet,  rocket,  or  compressor  air  exhaust  about  the 
axis  desired. 

Resolvent  matrix:  denoted  by  <J>(s);  the  matrix  which, 
when  multiplied  by  the  control  vector  and  the  Laplace  of 
the  input,  yields  the  Laplace  of  the  output  of  the  system. 

Root  locus:  a method  of  graphing  the  roots  of  a 
characteristic  equation  as  a function  of  gain;  useful  in 
stability  analysis. 
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Stability,  dynamic!  the  quality  of  a system  that 
determines  whether  its  response  to  any  input  will  be 
bounded  or  unbounded. 

Stability,  static:  the  tendency  of  a system  following 
displacement  from  a rest  position,  considered  oositive  if 
the  initial  tendency  is  to  return  to  the  original  position. 

State  variable:  a variable  which  determines  the  state 
of  the  system?  the  state  variables  of  a dynamic  system  are 
the  smallest  set  of  variables  which  completely  determines 
the  behavior  of  the  system  for  any  time  t £ 0 . 

Steady  s^.ate  response:  that  pare  of  the  total  response 
which  does  not  approach  zero  as  time  accroaches  infinity. 

Terminal  area  energy  management:  a dynamic  program 
using  inputs  of  orbiter  position  and  velocity  to  guide 
the  vehicle  to  its  final  approach  path,  minimizing  the 
heat  load  when  possible. 

Time  constant:  inverse  of  frequency. 

Time-invariant  differential  equation:  a differential 
equation  inwhich  none  of  the  terms  depends  explicitly 
upon  the  independent  variable,  t.^ 

Time  response:  a measure  of  the  output  of  the  system 
as  a function  of  time. 

3 Only  linear,  time-invariant  differential  equations  will 
be  considered  in  the  development  of  the  equations  of 
motion  in  this  study. 
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Transfer  function:  the  ratio  of  the  Laplace  transform 
of  the  output  (response)  of  a system  to  the  Laplace 
transform  of  the  input  (farcins'  function)  under  the 
assumption  that  the  initial  conditions  are  zero. 

Transfer  matrix:  the  matrix  of  transfer  functions 
for  a multiple  input/output  system. 

Transient  response:  that  part  of  the  total  response 
which  approaches  zero  as  time  approaches  infinity. 

Undamped  natural  frequency:  the  frequency  at  which 
a system  would  vibrate  were  it  not  damped. 

Unit  impulse:  an  input  whose  masnitude  approaches 
infinity  and  whose  duration  approaches  zero,  but  whose 
product  of  masnitude  times  duration  equals  one. 


A slossary  of  symbols  follows: 
a lift  curve  slope, 

dal 

b wing  span 

c mean  aerodynamic  chord 

Ct  * variation  in  lift  with  resoect  to  angle  of 

attack;  Generally  linear  at  angles  of  attack  below 
the  stall  region;  coefficient  of 
Cy  , variation  of  Ditching  moment  about  the  center 

* 2*x 

of  gravity  with  respect  to  angle  of  attack;  must  be 
negative  for  positive  static  stability;  coefficient  of 
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Zy0  pitching  moment  that  remains  when  lift  forces 
equal  zero,  coefficient  of 

e spar,  efficiency  factor,  relative  to  the  elliptical 
wina  and  lift  distribution 

f any  function,  time  domain,  as  f(t) 

? any  function,  Laplace  domain,  as  F(s) 

2 transfer  function,  time  domain,  as  s(t) 

transfer  function,  Laplace  domain,  as  3(s) 

I moment  of  inertia,  subscripted  as  to  axis  about  which 
it  is  taken 
K 2a  in 

<jC  Laplace  transform  of,  as  o££f(t)3 

cC  ' inverse  Laplace  transform  of,  as  J”'[F(s)l 

L,v,t  moments  about  X-,  Y-,  and  Z-axes,  respectively 

p,q,r  perturbations  in  roll,  pitch,  and  yaw  rates 

p,q,r  perturbations  in  roll,  pitch,  and  yaw  accelerations 

P,3,R  rates  of  roll,  pitch,  and  yaw 

S win2  area 

u,v,w  perturbations  in  forward,  side,  and  vertical 
velocities 

u,v,w  perturbations  in  forward,  side,  and  vertical 
accelerations 

U,V,W  velocities  in  forward,  side,  and  vertical  directions 
x,y,z  distances  al one:  X-,  Y-,  and  Z-axes 
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X,Y,Z  axis  system,  right-handed ; also  force  components 
in  X,  Y,  and  Z directions 
C<  angle  of  attack 

<5  Dirac  delta,  unit  impulse  function  cT(t) 

St  deflection,  elevon 
f damping  ratio 

Q'  § pitch  angle,  rate  perturbations 

0 pitch  angle,  relative  to  earth,  positive  nose  up 
f>  air  density 

<£  roll  angle,  relative  to  earth,  positive  right  wins- 
down  (see  also  resolvent  matrix) 
u>  frequency 

In  order  to  derive  the  equations  of  motion  of  a 
system,  a reference  frame  must  be  decided  upon.  The 
most  natural  of  these  would  be  a frame  in  which  the 
axes  were  aligned  with  the  physical  structure  of  the 
vehicle,  one  which  is  called  the  body-fixed  axis  system. 
3uch  a system  is  depicted  in  Figure  2-1 . Components  of 
velocities,  forces  and  moments,  and  distances  relative 
to  their  respective  axial  motions  are  listed  in  Table 
2-1 . 

In  the  development  of  the  equations  of  motion  for  a 
vehicle,  this  reference  frame  has  several  advantages. 

The  rotary  inertial  properties  are  constant  (assuming 
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motion 


velocity 


applied  forces 
and  moments  distance 
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constant  mass)  so  that  any  derivatives  with  respect  to 
time  are  zero.  The  pilot  senses  forces  and  moments  with 
respect  to  the  body-fixed  frame,  calculates  his  reactions 
to  vehicular  motion  and  to  his  changing  external  environ- 
ment in  terms  of  the  body-fixed  axes,  and  generally 
measures  his  position  (attitude),  velocities,  and  accel- 
erations from  a body-fixed  system.  This  is  especially 
true  when  he  is  dealing  with  primarily  relative  motion, 
as  in  air  combat  maneuvering. 


An  al 
space-fixed 
relative  to 
flat  earth, 


ternative  to  the  body-fixed  syst 
system  where  forces  and  moments 
an  axis  system  oriented,  for  ex 
as  in  Figure  2-2. 
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This  system,  however,  has  the  very  undesirable 
feature  of  requiring  inclusion  in  the  equations  of 
motion  the  variations  of  moments  of  inertia  of  a 
vehicle  as  the  vehicle  changes  its  velocity  vector 
and/or  attitude.  This  system  would  be  ideal  for  an 
aircraft  which  changed  neither,  but  such  is  r.ot  the 
case  in  this  study. 

Corrections  must  also  be  made  to  the  equations  of 
motion  in  the  body-fixed  axis  system,  but  these  corrections 
are  limited  to  inertial  forces  and  accelerations  and  do 
not  involve  changing  the  apparent  physical  characteristics 
o^  the  vehicle.  Further  simplifications  are  enabled  by 
this  choice  of  systems;  one,  the  small  ane-le  assumption 
is  made  possible  leading  to  the  linearization  of  the 
equations  of  motion  (see  Chapter  V);  two,  the  symmetry 
of  the  aircraft  across  its  X-Z  plane  eliminates  several 
terms  in  the  equations. 

The  body-fixed  axis  system  is  clearly  the  better 
selection  for  the  development  of  the  equations  of  motion 
of  the  vehicle. 

A definition  of  "flight  control"  is  necessary  in 
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order  to  identify  those  parameters  of  the  vehicle  with 
which  the  study  is  concerned.  The  definition  is  bor- 
rowed from  r/!cRuer,  Ashkenas,  and  Graham5  » 


"Control!  the  development  and  application  to  a 

vehicle  of  appropriate  forces  and  moments 
that 

1 . establish  some  equilibrium  state  of 
vehicle  motion  (operating  point  control) 

2.  restore  a disturbed  vehicle  to  its 
equilibrium  (operating  point)  state 
and/or  regulate,  within  desired  limits, 
its  departure  from  operating  point 
conditions  (stabilization)." 


This  definition  includes  the  concept  of  stability  and 
effectively  excludes  the  domain  of  guidance. 

In  light  aircraft,  all  functions  of  control  (turning-, 
climbing,  accelerating)  are  performed  by  the  pilot  using 
the  control  wheel  (or  stick),  the  throttle,  and  various 
other  aerodynamic  controls,  usually  the  flap  handle,  which 
are  directly  connected  to  the  device  which  they  operate. 

For  example,  the  pilot  desires  to  pull  up;  he  uses  some 
degree  of  back  stick,  the  elevator  responds  immediately 
and  in  proportion  to  the  amount  of  stick  travel,  and  a 
pitch  rate  is  generated.  The  pilot's  environment  remains 


5 Ibid. 
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relatively  constant  throughout  the  flight  envelope  of  the 
airplane;  a ten  pound  pull  on  the  stick  will  give,  for 
example,  1.3  g's  at  a given  airspeed  in  every  case. 

Even  for  higher  performance  aircraft  operating  in 
the  dynamic  environment  of  air-to-air  combat,  this  is  an 
entirely  satisfactory  system.^  Even  though  a given  ele- 
vator deflection  will  not  create  the  same  pitch  rate  at 
20000  feet  as  it  does  at  sea  level,  the  pilot  quickly 
compensates  for  this  by  modifying  control  stick  position. 
However,  maintaining  a particular  pitch  angle  can  better 
be  oerformed  by  an  automatic  control  system.  Using  a set 
of  accelerometers  to  sense  minute  deviations  from  an 
operating  (trim)  condition,  an  automatic  control  system 
can  command  the  rapid  and  precise  control  deflections 
necessary  to  maintain  the  vehicle  near  the  desired  flight 
condition.  A pilot  attempting  to  maintain  a given  pitch 
angle  will  find  his  sensory  threshold  much  higher  than 
(and  his  response  on  the  controls  much  less  precise  than) 
an  automatic  system.  This  type  of  system  is  installed  on 
many  aircraft  as  an  auto-pilot,  engaged  by  the  pilot  when 
desired,  to  maintain  certain  flight  conditions  such  as 

^ Present-day  fighters  have  flight  control  systems  aug- 
mented by  2-  or  3-axis  stabilization  to  minimize  gust 
response  and  other  oscillations. 
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heading:  and  altitude 


The  development  and  utilization  of  control  systems 
may  be  taken  one  step  further.  When  a vehicle,  such  as 
the  orbiter,  is  required  to  fly  in  vastly  different 
environments  in  the  course  of  a mission,  control  stick 
inputs  to  obtain  a given  maneuver  may  require  very  dif- 
ferent 'control  deflections.  Variations  in  dynamic  pres- 
sure, Mach  number,  and  altitude  create  changing  values 
of  control  effectiveness,  locations  of  aerodynamic  forces, 
even  stability.  The  automatic  control  system  can  use 
inputs  such  as  air  density,  dynamic  pressure,  Mach  number, 
and  existing  control  deflection  to  compute  precisely  how 
much  deflection  is  required  to  perform  a particular  maneu- 
ver and  how  rapidly  the  deflection  should  be  applied. 

The  system  may  also  consider  variations  in  stability  under 
different  flight  conditions.  It  will  then  command  that 
deflection,  measure  the  response,  and  make  further  cor- 
rections if  necessary,  with  far  greater  accuracy  and 
efficiency  than  a human  pilot. 

This  then  is  the  principle  by  and  for  which  the 
fly-by-wire  system  was  developed.  In  an  environment  of 
changing  flight  conditions  and  stability  where  high 
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control  precision  is  required,  the  automatic  system 
sigr.if icantly  outperforms  the  best  of  pilots. 
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III.  MATHEMATICAL  YODELS  OF  LINEAR  SYSTEM  ELEMENTS 


T’he  Laplace  trar.sforrr  is  a major  tool  in  the  study 
and  analysis  of  automatic  control  systems.  It  is  the 
procedure  by  which  time  domain  performance  is  connected 
to  frequency  response  and  results  from  the  correspondence 
betwee-  the  transfer  ^u^ctio’'-'  a^d  the  transient  resoor.se. 
The  transfer  function  itself  is  the  ratio  of  the  Laplace 
tra-.s^orn  or  a system's  output  to  the  transform  of  its 
input,  assuming-  that  all  initial  conditions  are  zero. 

A graphic  representation  of  such  an  input-output  rela- 
tionship is  provided  by  the  block  diagram.  In  the  block 
diagram,  of  which  the  transfer  fur.ctior.(s)  is  the  heart, 
inputs  and  feedbacks  are  represented  by  their  Laclace 
transforms,  and  the  paths  which  the  signals  follow  are 
depicted. 


Several  simplifications  in  computation  occur  when 
using  the  Laplace  transform,  among  which  are  ( 1 ) differ- 
entiation in  the  time  domain  is  represented  by  multipli- 
cation by  s in  the  Laclace  (frequency)  domain;  (2)  inte- 
gration becomes  division  by  s (both  with  regard  for  initial 
conditions);  and  (3)  convolution,  or  complex  multiplication, 
becomes  simple  multiplication. 


The  equation  linking  the  time  domain  to  * he  frequency 
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domain  is 

tv  c* 

lU<t)\  4 Fl»)  = & /e-’1  (It)  At  • ft'1  fit)  it 

t,  0* 

When  the  equation  of  motion  for  a system  can  be 
expressed  as  a linear  differential  equation,  the 
Laolace  transform  reduces  the  differential  equation 
to  an  algebraic  equation  in  the  Laplace  variable,  s. 
To  illustrate,  a damped  spring-mass  system  will  be 
used . 


/•  / / ^ ^ / / / / / / ^ ^ / < / / / 


x(t) 


Pismire  3-1 
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Here,  my  + cy  + ky  = x(t).  Each  term  is  replaced  by  its 
Laplace  transform: 

- « (»"«»)  - ,,i.)  - *j{»  ] 

tj  — c - iwl 

tj  — vtiO 
M)  -*  X(») 

When  coefficients  of  like  powers  of  s are  collected,  the 
equation  becomes 

+ cvjCo)  - XCs) 

The  appropriate  values  of  the  initial  conditions  are 
entered,  and  a simple  quadratic  equation  in  s has  replaced 
a second  order  differential  equation. 


Recalling  the  definition  of  a transfer  function,  all 
initial  conditions  are  set  equal  to  zero,  and  the  ratio 
Y(s)/X(s)  is  found : 

[rv^Cs)]  s"  * = X(0 

[rvvs^  c*  -v  W-l'fCs)  = XCO 


YCO i 

XCO  + cs  * k 


GCO  (3-2) 


where  G(s)  is  a common  symbol  for  the  transfer  function. 
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Figure  3~2  is  the  block  diagram  of  this  system  where 
X(s)  is  the  input  or  forcing  function,  Y(s)  is  the  output, 
and  G(s)  is  the  transfer  function  (also  called  the  plant 
matrix  in  multiple  input/output  systems).  In  this  case, 
G(s)  has  been  derived  from  mass  properties  of  the  system, 
the  degree  of  damping,  the  spring  stiffness,  and  funda- 
mental laws  of  motion. 


X(  s ) 


Figure  3-2 

In  modeling  a system,  G(s)  will  be  found  in  a manner 
analogous  to  the  above.  All  that  remains  i:  order  oo 
find  the  response  to  any  input  (neglecting  feedback  for 
the  moment)  is  to  multiply  the  transfer  function,  by  the 
Laplace  transform  of  the  input.  The  response  may  ce 
left  in  the  frequency  domain,  which  is  the  usual  case  in 
stability  analysis,  or  transformed  back  into  the  time 
domain  to  study  the  dynamic  characteristics  of  the 
system,  e.g,  rise  and  settling  times,  maximum  overshoot,  etc. 

Some  of  the  more  common  inputs  are  listed  in  Table 
3-1  with  their  time-  and  frequency-domain  representations . 
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The  transformations  are  linear:  that  is,  if  F(s)  is  the 
transform  of  f(t),  then  a»F(s)  is  the  transform  of  a*f(t)0 


time 
x(  t) 

frequency 

X(s) 

graph 

unit  impulse 

6 it) 

1 

J 

• 

unit  doublet 

<51 1) 

s 

J 

— 

unit  step 

J_ 

s 

J 

ramp 

t 

< 

s' 

y 

sinusoidal 

S\n 

_ \ v 

S -4- 

Table  3-1 

It  can  be  seen  in  the  first  column  that  the  unit 
step  is  the  time  derivative  of  the  ramp  function,  and  the 
unit  impulse  (represented  by  the  Dirac  delta)  is  the  deri- 
vative of  the  unit  step.  If  the  latter  is  not  obvious, 
consider  that  the  derivative  of  a function  is  the  slope 
of  that  function  at  every  point  where  the  function  is  con- 
tinuous, At  t=0,  the  slope  of  the  unit  step  approaches 


32 


infinity,  as  does  the  magnitude  of  the  unit  impulse.  At 
t=0+,  the  slope  of  the  unit  step  is  zero,  as  is  the  magni- 
tude of  the  impulse  function. 


Since  an  impulse  function  is  an  idealization,  it 
must  be  approximated  in  the  real  world,  When  finding  the 
Laplace  transform  of  this  approximation,  any  arbitrary 
magnitude  may  be  assigned  to  the  "impulse”,  as  long  as 
its  duration  is  small  compared  to  the  time  constant  of 
the  system.  If  the  duration  of  the  impulse  is  taken  as 
the  inverse  of  the  magnitude,  the  area  under  its  curve 
will  equal  one,  her.ce  the  unit  impulse. 


Unless  otherwise  depicted,  successive  terms  along 
a block  diagram  path  are  multiplied  together;  witness  the 
correlation  between  equation  (3-2)  and  Figure  3-2.  This 
is  true  whether  dealing  with  inputs,  successive  transfer 
functions,  or  other  operations  such  as  integration,  delay, 
or  differentiation.  Diagramatically, 


X(s) 


( s ) 


Go  ( s ) 


Y(s) 


is  identical  to 


X(  s )• 


•Y(s) 


The  resulting  transfer  function  G(s)  = G|  ( s ) « G2(  s ) , or 
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the  output  Y(s)  = X(  s ) • G-,  ( s ) • G2  ( s ) . 

A feedback  signal,  on  the  other  hand,  is  not  mul- 
tiplied but  rather  employs  a summing  junction  where 
certain  signals  downstream  in  the  diagram  are  fed  back 
and  added  to  (or  subtracted  from)  the  input  signal.  A 
typical  system  with  one  feedback  loop  is  depicted  in 
Figure  3-3«  —he  derivation  of  the  complete,  or  closed 
loop,  transfer  function  follows,. 


X(s)- 


^ \ ^ / 

G(s) 

CO 

H(s) 

. 

-»■  Y(  s ) 


Figure  3-3 


E(s)  represents  an  error  signal,  the  difference 
between  the  input  signal  and  the  feedback  signal.  3(s) 
is  called  the  open  loop  transfer  function?  it  is  the 
input  signal  after  being  modified  by  any  transfer  func- 
tions and  feedback  terms  up  to  the  summing  junction. 


In  this  case 

B(s)  = G( s ) »H( s ) • E( s ) . 
Since  Y(s)  = G(s)»E(s), 
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and  E(s)  = X(s)  - B(s)  = X(s)  - G( s ) • H( s ) • E( s ) 

then  E(s)  + G( s ) »H( s ) • E( s ) = X(s) 

or  E(  s ) jj.  + G(s)»H(s)j  = X(s) 

t?(s)  = X(  s ) 

1 + G(s)oH(s) 

Substituting  the  error  signal  above  into  Y(s)  = G(s)»E(s), 

the  transfer  function  of  the  system  due  to  feedback  be- 

comes  G.(s)  = Y(s)  _ G(s) 

L {SJ  X(s)  ‘ 1 + G(  s ) *H( s ) 0 


The  second  order  system  defined  by  the  damped  spring- 

mass  system  has  the  unique  definition  among  its  properties 

that  c K.  ■»_ 

^ and  7^  - where  J - damping  ratio 


and  = undamped  (natural)  frequency.  If  the  equation  of 
motion  is  divided  through  by  m,  the  result  is 


or, 


~ c.  • 

b * 7^ b ■* 

w 

/V\  b 

x(t) 

✓v\ 

? * 

= u>, 

Ct) 


where,  for  convenience,  the  input  function  x(t)  has  ab- 
sorbed a constant  l/k.  Talcing  the  Laplace  transform  as 
before  and  letting  the  initial  conditions  equal  zero, 
V(sT 


S*+  2 Twy.  5 -T  UJ„ 


t = XM 


V(Q  _ 

X (0  S * 2.  S * •-***•  . 
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If,  then,  the  transfer  function  of  a system  may  be  written 
in  the  format  b/(s^  + as  + b),  the  physical  characteristics 
of  the  system  are  immediately  evident. 

As  will  be  seen  in  the  analysis  of  the  orbiter,  the 
system  transfer  function  may  be  factored  into  two  such 
second  order  terms  yielding  at  once  the  periods  and  degrees 
of  damping  of  both  the  phugoid  and  the  short  period  motions. 
The  damping  term  offers  a further  clue  as  to  basic  or  in- 
herent stability  of  the  system.  When  the  roots  of  the 

quadratics  are  obtained  from  . r~z 

s=  * 5 J y -1  ^ the  value 

of  \ determines  whether  or  not  those  roots  lie  in  the  right 
half  plane  (RHP)  according  to  Table  3-2. 


value  of  J 

roots  are 

i<  o 

in  RHP 

1-0 

pure  imaginary 

o<  1 

complex  with  negative  real  parts 

r=  f 

equal,  negative,  and  real 

T > 1 

negative  and  real 

Table  3-2 

The  significance  of  a root  or  roots  in  the  right  half 
plane  will  be  discussed  in  the  section  on  stability. 
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IV.  STABILITY 


One  of  the  prime  indices  of  a vehicle's  perfor- 
mance is  the  measure  of  its  stability,  or  how,  how  much, 
and  how  quickly  it  responds  to  external  forces  or  control 
inputs . 

Two  types  of  stability  are  generally  defined,  static 
and  dynamic.  Static  stability  Is  the  vehicle’s  tendency 
to  return  to  its  original  flight  condition  after  being 
disturbed  from  that  condition.  A common  representation 
of  the  three  degrees  of  static  stability,  positive, 
negative,  and  neutral,  is  shown  in  Figure  4-1.  The  ball 
oossessing  positive  static  stability  will  return  to  its 
original  position  when  displaced,  the  ball  with  negative 
static  stability  will  continue  to  move  away  from  its 
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original  position,  and  the  ball  with  neutral  stability 
will  simply  remain  in  its  new  position. 

The  parallel  may  be  drawn  with  respect  to  an 
aircraft's  longitudinal  (pitch)  response  when  displaced 
from  trimmed  flight?  by  a control  incut  or  external 
force.  Figure  4-2  illustrates  the  longitudinal  forces 
acting  on  a vehicle  where  L is  the  sum  of  lift  forces 
produced  by  the  wing  and  tail  and  is  located  on  the 
aircraft  so  as  to  create  the  same  total  moment  about  the 
center  of  mass  as  would  the  individual  forces.  Cmo  is 
the  pitch  moment  remaining  when  the  lift  forces  are  zero. 


1. 


Figure  4-2 


? Trimmed  flight  is  defined  as  that  flight  condition 
where  the  sum  of  all  moments  acting  about  the  lon- 
gitudinal axis  is  zero,  i.e.,  = 0. 
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In  equilibrium  flight,  the  moment  produced  by  L 
is  balanced  by  Cm0.  and  no  pitch  rate  is  generated. 
Suppose  a srust  causes  a momentary  increase  in  the  anarle 
of  attack.  If  the  lift  curve  slope  is  positive  (that 
is,  a higher  anele  of  attack  generates  more  lift),  as 
are  all  modern  conventional  airfoils,  the  increased 
lift  will  cause  the  total  moment  about  the  center  of 
mass  to  be  nose  down,  or  stabilizing.  The  vehicle  will 
tend  to  return  to  its  original  condition. 

If  the  center  of  mass  should  move  aft  of  the  lift 
vector,  the  aircraft  could  still  be  stable  if  CmQ  were 
negative  (nose  down).  However,  if  that  sust  were  to 
hit  this  aircraft,  the  increased  lift  would  cause  a 
nose  UF  moment,  increasing  the  ane-le  of  attack,  further 
increasing  the  lift,  etc.  This  vehicle  then  would  be 
clearly  unstable. 

Three  quantities  are  necessary  to  this  discussion, 
the  lift  curve  slope,  the  pitching  moment  coefficient, 
examples  of  which  are  shown  in  Figure  4-3 , and  Cmo.  A 
vehicle  must  have  a positive  lift  curve  slope  ( jj-f  ),  a 
negative  pitching  moment  slope  (which  means  an  angle  of 
attack  above  the  trim  point  causes  a pitch  down  and  vice 
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versa),  and  a positive  Cmo  i-n  order  to  be  statically 
stable.  Cmw_  is  negative  when  the  aerodynamic  center 
of  the  vehicle  lies  behind  the  center  of  mass,  and  a 
positive  Cmo  is  designed  into  the  airframe. 

Given  that  a vehicle  has  positive  static  stabi- 
lity, its  dynamic  stability  may  be  studied'.  The 
ball  in  Figure  4-1 (a)  is  dynamically  stable;  if  dis- 
placed from  equilibrium,  it  would  return  to  the  bottom, 
probably  overshoot,  overshoot  again  coming-  back  (but  at 
a smaller  magnitude),  and  eventually  come  to  rest.  If 
there  were  no  friction,  the  system  would  be  dynamically 
neutrally  stable;  the  amount  of  overshoot  would  not  vary 
from  cycle  to  cycle.  If,  on  the  other  hand,  a magnetic 


3 An.  exception  to  this  requirement  for  positive  static 
stability  is  the  orbiter  itself.  Early  in  the  re- 
entry phase,  at  an  angle  of  attack  of  45°,  the  orciter 
is  laterally  unstable  due  to  the  blanking  of  the 
vertical  stabilizer.  However,  the  presence  of  a 
lightly  damped  Dutch  roll  gives  it  positive  dynamic 
stability. 
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device  were  installed  under  the  bowl  which  was  timed  to 
attract  the  ball  as  it  descended  and  repel  it  as  it  as- 
cended, enough  velocity  could  be  imparted  to  the  ball 
to  finally  shoot  it  over  the  side.  As  lone:  as  the  force 
supplied  by  the  magnet  was  by  any  degree  greater  than 
the  forces  of  friction,  the  ball  would  eventually  be  thrown 
from  the  bowl;  this  is  an  example  of  a statically  stable 
but  dynamically  unstable  system. 

Each  deerree  of  stability  can  be  exemplified  by  an 
aircraft.  A sharp  rap  on  the  control  stick  of  a trimmed 
aircraft  causes  a pitch  oscillation  which  quickly  damps 
out  and  disappears  (positive-stability) . Most  aircraft 
have  some  degree  of  Dutch  roll,  a continuous  but  mild 
combination  of  rolling  and  yawing  (neutral  stability). 

High  performance  aircraft,  whose  stability  may  be  mar- 
ginal under  certain  conditions,  can  be  susceptible  to 
FIO  (pilot  induced  oscillation)  where  the  pilot  plays 
the  role  of  the  aforementioned  magnetic  device  and 
amplifies  the  pitching  motion  of  the  vehicle,  with 
sometimes  catastrophic  results  (negative  dynamic  sta- 
bility). 

In  the  more  complicated  realm  of  automatic  control 
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systems  involving;  feedback,  high  precision  approach  and 
landing"  sequences,  etc.,  more  than  just  "positive"  or 
"negative"  as  regards  stability  must  be  known.  For  one 
reason  or  another,  it  may  be  desireable  to  utilize  a 
high  gain  for  a particular  feedback  signal,  but  there 
may  be  a limit  as  to  how  high  the  gain  may  be  before  it 
causes  the  system  to  become  unstable.  Certain  types  of 
incuts,  generally  frequency  related,  may  cause  large 
increases  in  the  magnitude  of  response  (resonance)  and 
may  require  a certain  amount  of  feedback  to  minimize 
these  excursions. 

Two  of  the  basic  tools  used  for  studying  the  effects 
of  feedback  on  a system  are  the  root  locus  graph  and  the 
Bode  plot. 

The  root  locus  is  a locus  of  roots;  the  roots  are 
those  of  the  characteristic  equation  of  the  system.  If 
the  real  parts  of  these  roots  are  all  negative,  then  the 
solution  of  the  differential  equation  describing  the 
system  will  involve  negative  (decaying)  exponents,  and 
the  system  will  be  stable.  This  is  the  rationale  be- 
hind requiring  roots  to  be  in  the  left  half  olar.e  for 
stability.  Feedback  gain  can  change  the  location  of 
these  roots  and  cause  a system  to  become  unstable. 
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An  example  is  in  order.  Consider  the  system  shown 
in  Figure  4~ 4 where  the  plant  transfer  function  is 

(3-2) 

(s+1 ) ( s+2 ) 

and  where  the  feedback  is  linear  and  variable  in  magni- 
tude. 


The  closed  loop  transfer  function  becomes 

( s-2 ) 

( s+1 )( s+2)  + K(s-2)  , 

and  the  characteristic  equation  is 

s2  + (3+K)s  + ( 2-2K ) = 0. 

The  roots,  for  varying  values  of  K,  are  given  in  Table  4-1, 


K 

root  #1 

root  #2 

0 

-2 

0.5 

-0.314- 

-3.186 

1.0 

0 

-4 

1.5 

0.212 

-4 . 71 2 

2.0 

0.372 

-5.372 

Table  4-1 


43 


Figure  4-5  is  a graph  of  the  root  locus  of  the  system. 

It  may  be  seen  that,  as  K exceeds  the  value  of  1,  one  of 
the  roots  moves  into  the  right  half  plane.  That  this 

t Jm, 


results  in  instability  is  shown  by  solving  the  corres- 
ponding differential  equation  for  a value  of  K > 1 . Let 
K=2;  the  characteristic  equation  becomes 

s2  + 5s  _ 2 = 0. 

Converting  this  equation  to  the  time  domain  by  taking 
the  inverse  Laplace  transform  yields 

y + 5y  - 2y  = 0, 

one  of  the  solutions  to  which  is  y = e®*3?2t#  Since 
y is  the  output,  it  may  be  seen  that  the  output  increases 
without  bound  as  the  time  increases.  The  system  is  un- 
stable for  this  value  of  feedback  gain. 


The  Bode  plot  depicts  how  the  system  will  respond 
to  a sinusoidal  input.  Gain  and  phase  margins,  which  are 
measures  of  relative  stability,  are  shown  directly  by  the 
plot,  as  is  the  magnitude  of  response  at  resonance.  By 
comparison  with  templates  of  the  proper  scale,  the  degree 
of  damping  may  be  read  from  the  graph. 


The  standard  representation  of  the  magnitude  of 
the  response  is  20*  log  |C-(  jco)|  where  jui  has  replaced  s 
in  the  transfer  function  expression.  This  magnitude  is 
plotted  against  the  frequency  of  the  input,  which  is  on 
a logarithmic  scale  so  that  wide  frequency  ranges  may 
be  plotted  with  accuracy.  The  phase  angle,  (S?  , is  like- 
wise plotted  against  the  logarithmically  scaled  frequen- 
cy. A stable,  time-invariant  system,  when  energized  by 
a sinusoidal  input,  will  have  a steady  state  output  of 
the  same  frequency  as  the  input  although  the  magnitude 
and  phase  angle  will  likely  be  different. 


Since  the  magnitude  of  the  response  is  plotted  in 
logarithmic  form,  finding  the  response  curve  of  a com- 
plicated transfer  function  such  as 


G(s)  = 


+ 3s  + 2 


s3  + 6s^  + 19s  + 13 
is  reduced  to  writing  the  function  as  the  product  of 


45 


first-order  factors,  finding  the  magnitude  of  each,  and 
then  adding;  (not  multiplying)  the  logarithms  of  the  mag- 
nitudes. From  this,  it  is  seen  that  an  increase  in  gain 
merely  raises  the  magnitude  curve  and  does  not  affect  the 
phase  angle. 


In  the  above  example,  writing  the  transfer  function 
as  a function  of  , 

-*•  ->-1  _ Q~U>V)  •» 


Letting  u>  = 0.1  yields 

U 


dv\  A 


l.oti. 
ivoi  «i 


o. 


Then  20*ln  | G(  jO . 1 ) I = -37.^38,  the  measure  of  magnitude  of 
the  response  when  the  input  frequency  is  0.1  radian  per 
second . 


The  phase  angle  is 


•w* 


CO 


l.SVl 
U ^ 


o 

O.'i'J.  . 


Computing  the  values  of  the  magnitude  and  phase 
angle  over  a wide  range  of  frequencies,  say  from  0.001 
to  10  radians  per  second,  enables  the  construction  of 
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he  3ode  diagram 


If  the  function  is  factored  into 


G(s) 


( s+1  ) ( s+2  ) 

( s+3+,12  ) ( s+3- j2  ) ( S+1  ) 


and  each  term  written  as  a function  of  jt*  , the  terms 

become  r r y' 

jti  + 1 jui+2  ^3  + j (0-2  )J  ^3  + j(uO+2)j  . 


The  magnitudes  au»  = 0,1  are,  respectively, 

1.005  2.002  0.282  0.273  0.995. 

The  logarithms  of  these  magnitudes  are 

0.005  0.69^  -1.266  -1.298  -0,005. 

which,  when  summed,  total  -1.870.  Multiplying  by  20,  as 
before,  yields  -37,^,  the  same  answer. 

The  sum  of  the  phase  angles  is 
5.71°  + 2.86°  - (-32.3°)  - 3^.99°  - 5.71°  = 0.17°, 
approximately  the  same  answer  as  before. 


The  factoring  method  is  desired  when  the  Bode 
plot  is  to  be  approximated  by  asymptotes,  or  straight 
lines  which  follow  the  trend  of  the  function.  It  does 
not  consider  damping.  The  former  method  is  entirely 
sufficient  for  plotting  the  exact  function  and  will  be 
the  one  used  in  this  study. 
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A third  type  of  graph  used  in  the  analysis  of  a 
system  is  the  time  response  graph,  or  time  history. 

Any  forcing  function  (not  just  sinusoidal)  for  which 
a Laplace  transform  exists  can  be  applied  to  a system, 
the  output  calculated,  and  the  output  inverted  into  the 
time  domain.  This  yields  a function  of  magnitude  with 
respect  to  time,  the  plot  of  whirh  depicts  period(s), 
overshoot( s ) , rise  times,  decay  times,  steady  state 
response,  and  gives  an  indication  of  the  degree  of 
damping  in  the  system.  Its  drawback  is  that  its  degree 
of  stability  is  r.ot  easily  derived,  although  it  is 
usually  apparent  whether  or  not  the  response  is  bounded. 

As  has  been  seen,  gain  will  move  the  locus  of 
roots  on  a root  locus  plot,  sometimes  changing  the 
stability  of  a system;  gain  will  move  the  Bode  magnitude 
plot  up  or  down,  changing  the  relative  stability.  When 
gain  is  applied  to  a system,  the  time  history  plot 
reflects  amplification  of  the  magnitude  of  the  response 
in  proportion  to  the  amount  of  gain,  i.e.,  twice  the  gain 
yields  twice  the  magnitude.9  For  example,  if  5°  of 
elevator  deflection  create  50  lbs  of  force  and  1000 

9 Since  we  are  dealing  with  linear  systems 
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ft-lbs  of  moment,  then  10°  of  deflection  will  generate 
100  lbs  of  force  and  2000  ft-lbs  of  moment?  the  pitch 
rate  created  by  the  latter  moment  will  be  double  that 
created  by  the  former.  However,  when  the  gain  to  be 
varied  is  in  a feedback  loop,  as  described  above  for 
the  root  locus  and  Bode  diagrams,  the  effects  on  the 
time  history  are  generally  unpredictable.  The  time 
history  is  therefore  of  limited  value  in  designing  a 
system;  its  asset  is  rather  in  analysis.  All  three 
types  of  plots  will  be  used  in  this  study. 
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V.  EQUATIONS  OF  MOTION 


In  confining  the  study  to  the  longitudinal  case, 
only  the  X and  Z forces  and  the  pitching  moment,  M, 
are  of  interest.  Each  of  these  may  have  derivatives  with 
respect  to  (1)  forward  velocity,  u,  (2)  pitch  angle,  q, 
(3)  vertical  velocity,  w,  (4)  vertical  acceleration,  w, 
and  (5)  elevon  deflection,  . These  derivatives  are 
labelled,  respectively,  Xu,  Xq,  Xw,  X^,  X^  , Zu,  etc. 

The  lower  case  subscripts  refer  to  perturbations  from 
steady  state. 

The  linearized  equations  of  motion  in  state  variable 
format,  in  terms  of  the  derivatives  defined  above,  are: 


Cl  = X.u  ♦‘X^LT  -V  - Xfcit 

K.  + ZjLT  ' (u.  *1y)v  * 

(.«->> 

I . «•  (M„ 

6 - * 
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The  state  variables  are  written  in  the  format 

u_ 


ur 


W = v 
6 


and  the  equations  of  motion  become 


-■fc 

U 

Z* 

(u.  * Z,v) 

0 

i. 

(h„  - Z.^Aw) 
c 

■*  7/w-  Ai Ir) 

[ j \ * (u  + zv)M.J 

i 

0 

0 

A 


Se 


J 


Or,  { *}  =•  [A]{*}  5e 


where  ^a"\  is  the  plant  matrix  and  is  called  the  control 
matrix . 


Some  justification  is  in  order  for  the  linearization 
of  the  system  equations  of  motion.  Assumptions  are  also 
made  to  simplify  the  development  and/or  to  define  the 
scope  of  the  study.  Reference  is  made  to  Figure  2-1  for 
sifin  convention  and  notation  in  a body-fixed  axis  systemi 
linear  and  angular  velocities  and  applied  forces  and  mo- 
ments are  identified. 
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The  first  assumption,  which  both  simplifies  the 
development  and  defines  the  environment,  is  that  the 
earth  is  fixed  in  space.  In  orbital  flight,  a universal 
inertial  frame  is  necessary  due  to  the  spinning  of  the 
earth  on  its  axis  and  to  the  more  pronounced  effects 
of  the  gravitational  fields  of  the  sun  and  moon  on  a 
spacecraft.  However,  for  the  relatively  short-term 
control  analysis  of  this  study,  a geo-inertial  frame 
of  reference  is  assumed. 

The  airframe  is  assumed  to  be  rigid.  This  enables 
the  description  of  vehicular  motion  as  a translation  of 
and  rotation  about  the  vehicle's  center  of  mass.  The 
method  of  application  of  the  thermal  protection  layer 
of  the  orbiter  requires  an  absolute  minimum  of  relative 
motion  between  and  within  various  sections  of  the  body, 
so  this  assumption  of  rigidity  is  especially  valid. 

The  mass  and  mass  distribution  of  the  vehicle  is 
assumed  to  be  constant.  Actually,  since  no  fuel  is 
being  burned  in  producing  thrust  and  no  stores  expended 
during  that  part  of  the  trajectory  under  study,  this, 
too,  is  an  especially  valid  assumption. 
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For  body-fixed  axes,  the  inertial  equations  of  motion 


(neglect 

ing 

eravit 

y)  are: 

= <w  0 + aw  - 

RV)  t 

X 

- (W  ■+  — 

= 

z> 

--  SI,  * (p-R")I..  = 

The  gravity 

force , 

acting  only  at  the 

vehicle ' 

s center 

of  mass, 

pro 

duces 

no  moments  but  does 

contribu 

ite  to  the 

external 

forces;  due  to  gravity, 

AX  : 

r 5<v(-  ©)  = 

w\c^  sin  0) 

h>Z  : 

: KWtj  Coi(-  6)  coi  - 

/Y\£^  Cj6S 

0 to:  4? 

The  corabinat 

ion  of 

the  inertial  and  gravitatio 

nal  force; 

on  the  vehic 

le  yie 

Ids : 

^Swv.  0^) 

X = 

<v\  (tj  -*•  SVJ  - RV  + 

z ~ 

(w  «■  pv  - a\J  - 

0)  tfiS  ' 

s) 

M.  -- 

- 

- ' 

pX. 

Dif feren' 

tiating  thi 

ese  equations  yield: 

3 equations  of 

perturbed  mo 

tion  (where,  for  example, 

dU  = u ) : 

dX  * 

> - w0 

4-  u”  -V.n  - 

' + C<^C03 

®.)e] 

a -- 

orv  \ 

w-  +•  V« 

■vT  -\J»  ^ - u- 

( 

4 

■ Oii  ©,  ) <P  + 

( ^ c*> &„)  g"] 

AM  - 

3 +• 

, * - at  P»  If  ) 
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where  the  quantities  subscripted  zero  are  the  steady 
state  (trimmed)  values.  To  arrive  at  these  equations, 

U was  replaced  by  U0  + u,  P by  P0  + p,  etc,  and  another 
assumption  was  madej  the  disturbances  from  steady  flight 
are  small  enough  so  that  the  sine  of  an  angle  may  be  re- 
placed by  the  angle  itself  and  the  cosine  of  an  angle 
equals  one.  Any  products  of  perturbations  are  neglected. 
The  "small  perturbation"  concept  has  linearized  the 
transcendental  equations  of  the  system. 

Because  of  the  symmetry  of  the  airframe  in  the  X-Z 
plane,  and  because  there  are  no  aerodynamic  asymmetries 
such  as  propeller  slipstream,  it  is  assumed  that  lon- 
gitudinal forces  and  moments  due  to  lateral  perturba- 
tions about  the  trim  condition  are  negligible.  Pitch 
rate  at  trim,  30,  is  defined  to  be  zero.  A further 
simplification,  which  now  has  no  effect  on  the  system, 
is  that  P0  = R0  = VQ  = <$>0  = zero.  This  yields 

ax  - ~ 

VL  - ^ - -U.  ^~[ 

dfA  = \ 

Since  the  aerodynamic  forces  and  moments  on  the 
vehicle  are  dependent  mainly  upon  the  velocity  of  the 
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vehicle  with  respect  to  the  air  and  not  upon  accelera- 
tion of  the  air  mass,  all  derivatives  with  respect  to 
acceleration  are  neglected.  The  exception  is  w,  an 
important  one  as  it  corresponds  to  the  time  lag  between 
a vertical  velocity  perturbation  at  the  leading  edge 
and  its  effect  at  the  trailing  edge.  It  would  not  exist 
if  downwash  did  not  exist,  but  by  this  very  fact,  its 
magnitude  cannot  be  disregarded.  The  flow  is  now 
assumed  to  be  quasi-steady. 

Altitude  perturbations,  being  assumed  small,  allow 
the  neglect  of  air  density  and  temperature  perturbations. 

If  the  body  axes  are  now  rotated  so  that  the  x-axis 
points  into  the  relative  wind,  W0  may  be  set  equal  to 
zero.  This  orientation  forms  a stability  axis  system. 
(Since  small  angle  assumptions  have  been  made,  pertur- 
bations can  still  be  measured  in  the  body-fixed  axis 
system. ) Applying  this  rotation  of  axes  to  equations 
(5-1)  and  taking  the  Laplace  transform  of  each  yields 

^ - £-r)  W 

* Mvr)  U-  = Mft<y€ 
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A vector  diagram  of  the  forces  acting  on  the  ve- 
hicle will  aid  in  the  discussion  of  the  derivative  terms 
in  the  above  equations;  see  Figure  5-1.  Note  that  there 
is  no  thrust  force  present. 


Figure  5-1 

The  sums  of  forces  in  the  x-  and  z-directions  is, 
for  small  <*  X = Urn  a - Du»  a = Lot  - D 

Z » - Ql  - Da) 

The  velocity,  V,  is  assumed  to  be  equal  to  U0. 
This  assumption  is  justified  by  the  following  treat- 


56 


Since  u,  v,  and  w are  perturbations  and  are  small  v/hen 
compared  to  U0,  their  products  may  be  neglected! 


= \J V o2  + 2U0u  = U0  yj  1 + u/lJ0  2 l 


U0 


Since  £0t=  , then 


dX  / o,  , « zt>\  _ _ u 


where 

<Xe  is 
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to  be  zero.  Us 

X--- 

♦n  "51*. 
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X.  --  ~ 

t-  (tjx-*Sp  + 

5m  v.  act  + 

pSU.  /tr.  2>C, 

” — 1 

m V 
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. „ / au  . al>  k\  2>l 

^ r Su  * = ^ **  * ^ *)  = " ^ 
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, then 
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«v\  V 

’ *£• 

where 

the  moment  = M = f> 

Mu  * 

pUSc  / 
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Using  the  same  substitution  and  differentiation 
techniques  as  above,  the  rest  of  the  force  derivatives 
are  found. 


The  q-derivatives i 


\T 

fU.S£ 

sCp 

4 AO 

c)  (%£/acr.) 

^ ' 

fU.Sc 

*S><Lu 

4m 

^ 0»Vxo.) 

\ - 

■>  OVw.) 

w-deri 

vatives : 

c>Cp  'N 
~ b*  / 

l» 

* 

t 

( 

C>Ct-  f \ 

Si  * 

Mvj-  - 

pV5,Sc 

^C,v, 

"hot- 

w-derivatives : 

x*  - 

S Cp 

xSc 

Zw 

K^u.) 

il  - 

*iT.  _ 

10 


e>("%Vj0) 


10The  r.on-dimensior.alization  procedure  is  discussed  later 
in  this  chapter. 
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The  ^-derivatives! 


— 

X*  = 

^0  o 

Is.  -- 

_ 

^ - 

ptJ.^Sc 

at. 

Since  the  forces  and  moments  are  proportional  to, 
respectively,  mass  and  moment  of  inertia,  the  stability 
derivatives  are  non-dimensionalized  by  dividing:  the 
force/moment  equations  by  these  quantities;  hence,  the 
appearance  of  l/m  and  l/ly  in  the  results. 

Some  of  the  derivatives  have  particular  signifi- 
cance or  are  named  accordin.gr  to  their  effect  on  the 
motion  of  the  system.  Others  may  be  neglected  for 
various  reasons,  A discussion  of  these  items  is  in 
order. 

The  quantity  is  called  the  speed  damping  deri- 

vative. 

e>tu  . 

The  quantity  is  generally  small  except  at 
transonic  speeds,  the  subject  of  this  paper. 

The  quantity  results  from  aeroelastic  effects, 
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which  are  being  neglected,  and  from  compressibility 
effects,  which  are  not.  It,  too,  is  largest  at 
transonic  speeds. 

Cmq  is  damping  in  pitch  and  is  negative  for  a 
stable  vehicle.  Its  primary  effect  is  on  short  period 
damping. 

The  quantity  ^ is  the  static  stability  deriva- 
tive, the  prime  indicator  of  static  longitudinal  sta- 
bility. liven  an  increase  in  ar.gle-of-attack , the 
increased  lift  on  the  wing  (normally,  a tail  contri- 
bution is  dominant)  will  cause  a negative  pitching 
moment  (nose  down)  about  the  center  of  mass  of  the 
vehicle.  If  this  is  the  case,  the  center  of  mass  must 
be  forward  of  the  aerodynamic  center  of  the  wing,  and 
the  vehicle  is  statically  stable.  Considering  the 
tailless  orbiter,  it  is  expected  that  Cm  will  be  small 
compared  to  a conventional  aircraft,  becoming  somewhat 
more  negative  (more  stable)  as  the  vehicle  becomes 
supersonic  and  the  aerodynamic  center  moves  further 
aft  on  the  wing. 

The  quantity  ^ is  the  lift  curve  slope  and  is 
positive  and  generally  linear  at  angles-of-attack 
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below  stall 


Since  a change  in  drag  on  a horizontal  tail  is  the 
main  contributor  to  a change  in  X„,  and  since  this  drag 
increment  is  small  anyway  compared  to  the  total  draff, 

Xw  is  neglected  in  further  discussion. 

Because  of  the  time  lag  between  a change  in  angle- 

of-attack  at  the  leading  edge  and  the  change  in  pressure 

*CU 

distribution  at  the  trailing  edge,  is  produced. 

However,  it  is  a negligible  quantity. 

The  quantity  will  increase  short  period  damping 

when  it  is  neffative,  as  is  normally  the  case.  It  is 

, . . , . , 

ron-dimensioralized  as  — 

K*Vzv.)  ' 

7>Ci.  ICt, 

The  quantities  and  are  both  significant  in 

tailless  vehicles,  the  first  being  always  positive,  the 
second  smaller  and  sometimes  negative.  Their  ratio  may 
be  likened  to  L/D,  Cmr^  is  elevator  control  effective- 
ness and,  as  the  name  implies,  is  quite  important.  It 
determines  how  much  pitching  moment  will  be  generated 
for  a given  elevon  deflection.  It  is  sensitive  to 
center  of  mass  travel,  aerodynamic  center  travel  (as 
when  a vehicle  becomes  supersonic),  and  CLmax*  It  is 
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a negative  quantity  for  elevons  aft  of  the  center  of 
rravity . 

Since  « = tan_1(w/u0)  = w/’J0,  the  terms  oc  and  w may 
be  interchanged  (as  may  their  derivatives)  with  due 
respect  for  the  constant  of  proportionality# 

The  non-dimensional  stability  derivatives  intro- 
duced on  pa<re  61  will  be  used  wherever  practical. 

In  accordance  with  7!cr?uer's  definitions,  they  are 
called  basic  derivatives  because  they  do  not  involve 
any  inertial  quantities.  They  are  called  non-dimensional 
because  they  involve  force  and  moment  coefficients 
rather  than  the  forces  and  moments  themselves.  They 
are  useful  in  that  study  of  the  vehicle  under  different 
flight  conditions  is  most  easily  correlated.  It  may  be 
seen  that,  for  example,  is  the  non-dimensional 

form  of  . 

The  other  type  of  derivative,  of  which  Xu  is  an 
example,  is  called  a dimensional  stability  derivative 
parameter.  It  will  be  used  in  the  study  of  airframe 
transfer  functions. 
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VI.  COMPUTER  METHODS  OF  ANALYSIS 


The  concept  of  state  variable  analysis  was  made 
possible  with  the  advent  and  refinement  of  the  nigh 
speed  digital  computer.  The  arithmetical  operations 
involved  in  dealing  with  just  a 4x4  matrix  are  copious; 
they  include  inverting,  factoring,  solving  the  .equations 
the  matrix  represents  (which  are  often  differential 
equations),  and  applying  the  results  to  state  variable 
procedures  for  obtaining  such  things  as  the  resolvent 
matrix,  the  state  transition  matrix,  and  the  time  res- 
ponse. In  addition,  most  of  these  results  must  be  com- 
puted for  many  instances  of  time,  say  every  0.1  second 
for  200  seconds,  or  for  a wide  range  of  frequencies  in 
order  that  meaningful  data  may  be  plotted  and  analyzed. 

3efore  the  computer  methods  can  be  presented,  a 
description  of  the  underlying  state  variable  theory  and 
mathematics  is  necessary. 

Conventional  control  theory,  exemplified  by  the 
damped  spring  mass  system  in  chapter  III  (pp  29  et 
seq„),  is  adequate  for  sinele-input-single-output 
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operations  but  is  wholly  inadequate  for  multiple-input- 


multiple-output  systems.  Each  additional  parameter 
introduced  as  an  input  or  output  increases  the  number 
of  mathematical  operations  many-fold. 

The  set  of  state  variables  must  first  be  specified. 
In  accordance  with  the  equations  of  motion  (p  53  ) , the 
set  will  be  (u,w,q,0).  The  state  vector  is  seen  to  be 


the  control  vector  is  u = , and  the  plant 


matrix  and  control  matrix  are  shown  as  [a] 
pectively  on  pas:e  5-2.  The  state  equation 
tern  is  then 


{xj  = \fo)  + {3}u  . 


and  ^ res- 
for  the  sys- 


The  Laplace  transform  is  used  to  solve  this  non- 
homoe'eneous  equation  of  statet 


s'Xlh  - *(<0  • [a3X<0  - £C(-s) 

-(MS')  * £<■»)  * §.w<) 

- £<“■>*  %'K,) 

y(0  - 
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Since  (si.- A)  = 5-  ^ A A- 

and  X*'  [(.sT-A)*']  5 1 + ^ + IT  ^ IT  + ’ " - 

then  XCO  - X [eAt]  x(6)  -r  X[cM]B15CO 

Taking  the  inverse  Laplace, 

*(t)  - eAtxC*i)  -*■  /e'^  ' fc»u(t)<lT  (6-1) 

~ o — 

where  T is  a dummy  variable. 

The  matrix  eAt  is  written  as  «$(t)  and  called  the 
statg  transition  matrix.  Its  Laplace  counterpart,  $(s), 
is  called  the  resolvent  matrix.  The  first  term  on  the 
right  hand  side  of  equation  (6-1)  is  the  time  response 
due  to  initial  conditions;  the  second  term  is  that  due 
to  the  input,  Bu. 


Just  as  in  the  single-input-single-output  system, 
where  Y(s)  = G(s)*X(s),  here,  X(s)  = G(s)*U(s),  but  G(s) 
is  now  a matrix  and  is  called  the  transfer  matrix.  The 
corresponding  closed  loop  system  is  depicted  in  Figure 


6-1. 


U(s) 

3(  s ' 


Go  ( s ) 

rm 

j 1 

lH(s)  k 

:>  X(s) 


Figure  6-1 
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As  was  developed  for  the  single  ir.put/output  case, 
B(s)  = H( s ) *X( s ) 

or,  dropping-  the  "(s)"  for  convenience, 


B = HG0E 

Since 

X = G0[u  - B] 

ii 

a 

o 

i 

n: 

then  X 

+ G0HX  = G0U 

t1  + 

G0h]x  = G0'J 

and 

x =[i  + g0hXg0u 

'Writing  X as  G’J  shows  the  closed  loop  transfer  matrix 
to  be  C(s)  = [i  + G0(s)*H(s)]  *G0(s).  1 

In  order  to  solve  the  equations  of  motion,  a 
computer  program  will  use  the  above  relationships  and 
some  form  of  numerical  integration.  The  programs 
written  by  IVelsa  and  Jones  were  used  in  this  study, 
the  main  two  of  which  were  BAS^AT  and  RTRESP. 


1 1 


Since,  in 
exercised 
that  both 
matrix  in 


general , [a][b]  * Ib]{a],  care  must  be 
in  developing  the  equations  by  insuring 
sides  are  pre-  or  post-multiplied  by  a 
the  same  order. 


J.L.  IVelsa  and  S.K.  Jones,  Computer  Programs  for 
Computational  Assistance  in  the  Study  of  Linear 
Control  Theory , 2nd  edition,  iVjcGraw-Hill , 1973 
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Given  a basic  matrix,  [a] , in  this  case  the  plant 
matrix,  BASMAT  produces 

a,  the  determinant  of  [a]  , det[_A] 

b.  the  inverse  of  La],  LaV' 

c.  the  characteristic  polynomial,  det(sI-A) 

d.  the  eigenvalues  of  [A},  X; 

e,  the  state  transition  matrix,  $(t)  = e*1 

• | 

f,  the  resolvent  matrix,  $(s)  = (sI-A) 

Once  the  state  transition  matrix  is  known,  the 
time  response  may  be  found  by  equation  (6-1)  whether  the 
input  is  in  the  form  of  an  initial  condition( s ) , a for- 
cing function,  or  both.  The  other  Melsa-Jones  program, 
RTRESP  (Rational  Time  RESPonse),  was  used  for  this 
purpose . 

After  the  mathematical  formulation  of  the  time 
response  is  obtained,  it  is  left  to  graph  the  results 
to  see  the  time  history  of  the  system.  In  order  to 
facilitate  this  operation,  a program  named  GRAPH  was 
written,  in  BASIC  language,  for  execution  on  the  HP9830 
computer.  See  Figure  6-2.  An  explanation  of  the  vari- 
ables and  operations  in  GRAPH  follows i 
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1 LI1’  r MJuRilii  H(  H 


l O Dll!  CC  ] ■.  r 1 1 Qi  9 j . RL  9 j 
20  FOR  1=1  TO  9 
30  CC  I 3=Pr  1 3=01  1 ]=Rt  i 3=0 
40  NEXT  I 

50  DISP  “HOW  MANY  REAL  ROOTS"? 

60  INPUT  F:1 

70  DISP  “HOW  MANY  COMPLEX  ROOTS"? 

SO  INPUT  R2 

SO  IF  R1=0  THEN  140 

100  FOR  1=1  TO  Rt 

110  DISP  "INPUT  REAL  ROOT  #"I 

120  INPUT  PC  I 3 

130  NEXT  I 

140  FOR  I =R 1+ 1 TO  R1+R2  STEP  2 

150  DISP  "REAL  I MAG  PARTS  OF  CPLX  ROOT  #“I? 

160  INPUT  OCI3»RCI3 
170  QC 1 + 1 J=C1 I 3 
ISO  RC 1 + 1 ] — -RT  I 3 
190  NEXT  I 
2O0  WRITE  <15»210> 

210  FORMAT  45“*" 

220  PRINT  "ROOTS  OF  CHAR  EQN  ARE" 

230  PRINT  " REAL  PART  I MAG  PART" 

240  FOR  1=1  TO  R1 
250  WRITE  <15»310>PC I 3»0 
260  NEXT  I 
270  PRINT 

280  FOR  1 =R 1 + 1 TO  R1+R2 
290  WRITE  <15»310>QC  I 3>RC  I 3 
300  NEXT  I 

310  FORMAT  4X  > F 1 0 . 5 » 4X  > F 1 0 . 5 

320  IF  R1=0  THEN  370 

330  FOR  1=1  TO  R1 

340  D I SP  " COEFF  OF  EXP  < " PC  13"  T ) " » 

350  INPUT  CC I 3 
360  NEXT  I 

370  FOR  I=R1+1  TO  R1+R2  STEP  2 
3S0  RC  I 3=ABSCRC  ID 

390  DISP  “COEFF  OF  EXP < M QC I ] ” T ) *CQS ( " RC I 3 " T > " ? 
4O0  INPUT  CC 1 3 

410  DISP  "COEFF  OF  EXP< "QC I ]"T>*SINt "RC I 3"T> " ? 
420  INPUT  CCI+13 
430  NEXT  I 

440  DISP  "MAXIMUM  MAGNITUDE  OF  RESPONSE"? 

450  INPUT  G9 


Fleur*  4*2 


M 


460  .SCALE  -5'  *..>-1.  1*G9>  l.  1 *G9 
•♦70  XRX 1 3 0.5.0*10 
4k:0  YAM  I S -0 . 8 * G9  4 i -G9  * G9 
490  LABEL  < *> 2. 5. 1 . 7. 0. 10- 7) 

500  FOR  X=0  TO  10  STEP  5 

5t0  plot  x.o.i 

520  CF'LOT  -2.5. -2 
530  LABEL  v770)X 
540  NEXT  X 

550  FOR  Y=-G9  TO  G9  STEP  2*G9 

560  PLOT  0 » Y » 1 

570  CF"  OT  -7*  -0.3 

580  LABEL  <780)Y 

590  NEXT  V 

600  FOR  T=0  TO  10  STEP  0.25 
610  Y=0 

* GO SUB  1080 

PLOT  T * Y 
t »0  NEXT  T 
650  PEN 

660  DISP  "RESET  PLOTTER  FOR  PHUGOID 
STOP 

6 SO  SCALE  -15.215. -1 . 1*G9> 1 . 1*G9 
6 90  XAXIS  0.25*0*200 
700  YAXIS  -8.G9/4.-G9.+G9 
710  LABEL  <*>2.5* 2.0.7- 10) 

720  FOR  X=0  TO  200  STEP  50 
730  PLOT  X.0.1 
740  CPLOT  -2. 5, -2  • 

?50  LABEL  <770)X 

760  NEXT  X 

7-0  FORMAT  F4.0 

*,'-0  FORMAT  F8.3 

'90  FOR  T=*0  TO  10  STEP  0.25 

800  Y=0 

810  GOSUB  1088 
820  PLOT  T.Y 
330  NEXT  T 

840  FOR  T=10  TO  50  STEP  1 
850  Y*0 

860  GOSUB  1080 
870  PLOT  T.Y 
880  NEXT  T 

890  FOR  T=50  TO  200  STEP  2 
900  Y*0 


Figure  6-2  (continued) 
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910  GO SUB  1030 
920  PLOT  T»Y 
930  NEXT  T 
940  PEN 

950  DISP  "INPUT  LETTER  SIZE  FOR  LABELLING" $ 

960  INPUT  Z 

970  LABEL  <*» Z» 1 . 7> O, 7/l0> 

980  DISP  "YOU  ARE  IN  THE  LETTER  MODE" 

990  LETTER 

I860  DISP  "ANY  MORE  LABELLING"? 

1010  INPUT  Z 

1020  IF  Z=1  THEN  950 

1030  DISP  "ANY  MORE  GRAPHS  WITH  THESE  ROOTS"? 

1040  INPUT  Z 

1050  IF  2*1  THEN  320 

I860  DISP  " THE  END " 

1070  STOP 

1080  REM  SUBROUTINE  TO  COMPUTE  MAGNITUDE  OF  RESPONSE 

1090  IF  R1*0  THEM  1140 

1100  FOR  1*1  TO  PI 

1110  IF  PC  I ]*T\ ' 5 THEN  1130 

1120  Y*Y+CC I ]*EV>  vPC I ]*T> 

1130  NEXT  I 

1140  FOR  I*R1+1  TO  R1+R2  STEP  2 

1150  IF  QC I ]*T  >225  THEN  1180 

1160  Y*Y+CC  1 ]*EXP<QC  I ]*T>*COS<RC  I ]*T> 

1170  Y*Y+CC 1 + 1 ]*EXP<QC I ]*T)*SIN<RC I 3*T) 

1130  NEXT  I 
1190  RETURN 
1200  END 


Plfm  6-2  (continued) 
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C^i  coefficients  of  the  exponential  and  sinusoidal 
terms,  the  sum  of  which  describe  the  time  res- 
ponse 

G9 ; maximum  absolute  magnitude  of  the  responses  used 
to  scale  and  label  the  time  history  plot 

I i counting  element 


P^;  purely  real  roots 

Qi s real  part  of  complex  roots 

Rj_«  imaginary  part  of  complex  roots 

R1 ; number  of  purely  real  roots 

R2 : number  of  complex  roots 

T;  time;  see  "X" 

X;  horizontal  (time)  axis 
Yi  vertical  (magnitude)  axis 

Z:  yes(l)/no(0)  value;  used  to  answer  various 

questions  which  the  program  asks  of  the  operator 


lines  1 0—9-0  s dimensions  C,  F,  Q,  and  R and  initializes 
their  value  at  zero 


11  50-3"!  o ; inputs  roots,  categorizes  them  as  real  or 
complex,  and  prints  results 

11  320-9-30;  inputs  coefficients  of  exponential  and 
sinusoidal  terms 

11  49-0-590;  scales,  draws,  and  labels  the  axes  of 
the  graph 

11  600-650;  plots  short  period  magnitude-vs-time 
graph,  using  subroutine  1080 

11  660-99-0;  plots  phugoid  graph,  using  subroutine  1080 
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11  950-1020:  allows  for  labeling; 

11  1030-1070i  allows  drawing  of  other  graphs  with 
the  same  roots 

11  1080-1190:  subroutine  which  computes  magnitude 
of  response  at  a given  time 

line  1200:  END  statement 

. The  values  of  the  coefficients  and  roots  (Cj_,  Pi,  Qj_, 
Ri)  are  obtained  from  the  program  RTRESP. 


The  text  Aircraft  Dynamics  and  Automatic  Control, 
referenced  previously,  contains  an  input/output  analy- 
sis of  an  early  operational  Air  Force  interceptor,  the 
F-89  Scorpion.  In  order  to  validate  the  GRAPH  prog-ram, 
this  example  was  run  through  the  3A3MAT  and  RTRESP  pro- 
grams and  the  results  plotted  using  GRAPH. 


From  the  given  stability  derivatives,  the  plant 

matrix  and  control  vector  were  constructed  in  accordance 

with  the  state  equations  on  page  50  5 for 

altitude  = 20000ft 
weight  = 305001b 
Mach  = 0.638 
TAS  = 660fps 


U = -0.0097 

Zu 

= -0.0955 

Mu 

= 0.0 

w = 0.0016 

Zw 

= -1.430 

Mw 

= -0.0235 

4c  = °*° 

z$t 

= -69.8 

Mw 

= -0.0013 

Mq 

= -1.920 

= -26.009 
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the  plant  matrix  £a]  was 


0.0097 
0.0955 
0.0001 24 
0.0 


0.0016 
-1 .430 
-0.0216 
0.0 


0.0 
660 , 0 
-2.778 
1 .0 


-32.2 

0.0 

0.0 

0.0 


and  the  control  vector  was 


0.0 

-69.8 

-26.009 

0.0 


The  input  to  the  system  consisted  of  an  elevator 
deflection  of  0.02  radian  commencing  at  t=0  with  a 
duration  of  one  second;  initial  conditions  were  equal 
to  zero.  The  text  results  are  shown  in  Figure  6-3. 


In  order  to  utilize  the  RTRESP  program,  a proper 
Laplace  transform  of  the  input  had  to  be  found.  The 
first  choice  was  a multiple  of  the  unit  impulse,  in 
this  case  0.02  <5  ( t ) . The  Laplace  of  this  input  is 
0.02(l/l).  However,  the  Melsa-Jones  program  requires 
that  the  denominator  of  the  Laplace  of  the  input  be 
one  or  more  orders  greater  that  that  of  the  numerator. 
Since  this  prerequisite  could  not  be  satisfied,  it  was 
decided  to  use  a step  function  input  and  take  the  deri- 
vative with  respect  to  time  of  the  output,  since  the 


73 


20  - 

t 


-20  - 


9 

(rod^ec) 


(ft/sec) 


-20  - i 


9 

(rcd/Vc) 


■I  t- 

ol-L 


I 


- 1 


I 


• • uw  » 

-1000  - -1000  - 


04- 

*. 

0 • 

(rod! 

■04  I- 


.04 


(rod) 


0 -• 


■.04 


6 io 

time , t { sec ) 

O ) Short  PerioJ 


-i.__I  i.-J.-L  1 i i I , . i.  I l . L - X _L 

o 50  100  150  200 

time  ,t  (sec) 

t)  PfiLgcid 


Figure  6-3 

impulse  is  the  derivative  of  the  step  function.  The 
results  of  the  system's  matrix  analysis  are  given  in 
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-ip-ure  6-4.  The  results  of  ;he  step  function  input  are 
show-,  in  Firure  6-5.  and  Table  6-1  shows  that  the  periods 
and  damping  ratios  of  the  short  period  and  the  phugoid 
agree  quite  closely  with  the  text  results. 


SHORT  PERIOD 

PHU00ID 

3A3IVAT 

text 

3 ASK AT 

text 

C0<r\ 

4.2693 

4.27 

0.0629 

0.0630 

T 

0.492^ 

0.493 

O 

• 

O 

0.0714 

Table  6-"1 


The  differentiation  of  the  results  of  the  step 
function  innut  yielded: 


let  a i = -2.104339  bx  = 3.714746 

a2  = -0.004510971  b2  = 0.06274506 


u = -0 , 599'48exp( ai t )cos( bi t ) 

-O.6781 76 exp ( ai t ) sin( bi t ) 
+0.599209exp(a2t)cos(b2t) 

+20. 06399 7exp( a2t)  sir.(b2t ) 


4 
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Figure  6-4  (continued) 


********************************************* 
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Pigure  6-4  (continued) 


RATIONAL  TIME  RESPONSE 

PROBLEM  IDENTIFICATION  - MCP.UER  EXAMPLE 
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Figure  6-5 


DENOMINATOR  POLYNOMIAL  OF  R<$)  - ASCEKOINS  PO'-.'FRS  OF 
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Figure  6-5  (continued) 


w = -1 . 336301  exp(  a-!  t )cos  ( bi it ) 

-92. 66498 5exp(ai t)sin(bi t) 
-0.005^06exp(a2t)cos(b2t) 

-0. 195852exp(a2t)sin(b2t) 

q = -0,520226exp(ai t)cos(bi t) 

+0.1 02383 exp (ai t)sin(bi t) 

+0. 00003 ^ exp ( a2t) cos ( b2t) 

+0 . 002468 exp ( a2t ) sin( b2t ) 

Q = +0.039993exp(ait)cos(bi it) 

-0.1 1 785b exp ( ai t) sin(bi t) 

-0. 0391 89exp(a2t) cos (b2bj 
-0 . 002065exp( a2t ) sin( 83 1 ) 

Figure  6-6  shows  the  time  history  of  this  system. 
The  magnitudes  of  vertical  velocity,  pitch  angle,  and 
to  a lesser  degree,  pitch  rate  were  decidedly  in  error. 
Furthermore,  the  initial  condition  of  zero  pitch  rate 
was  not  reflected  in  the  results. 

The  discrepancies  are  due  to  the  improper  descrip- 
tion of  the  input.  Clearly,  the  effect  of  a maximum 
elevator  deflection  acting  over  a short  period  of  time 
(A  t « 1)  is  different  from  a much  smaller  deflection 
for  one  second,  hence  the  magnitude  errors.  The 
initial  condition  error  is  due  to  the  indefinite 
integration  of  the  input.  When  the  input  is  integrated 
(multiplied  by  s),  the  denominator  loses  a root  at 
s=0  and  its  associated  eigenvalue.  The  time  response 
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Figure  6-6 


therefore  loses  an  exponential  term,  and  the  initial 
conditions  are  not  satisfied. 


The  undamped  natural  frequencies  and  damping  ratios 
are  calculated  from  the  roots  (eigenvalues)  of  the 
characteristic  equation  by  considering  the  roots  in 
vector  form. 


where  the  roots  are  expressed  as  X = a + jb  and  were 
calculated  by  BASMAT. 


The  second  attempt  to  duplicate  the  given  input 
function  was  to  evaluate  the  results  of  the  step  function 
(xO. 02)  at  t=l  and  use  these  values  as  initial  conditions 
for  a zero-input  system  commencing  at  t=l . The  values  of 
the  state  variables  at  t=l , calculated  from  the  RTRESP 
program  (see  Figure  6-5)  were 


84 


0.9938) 
j -21.8208  V 
\ -0.0035 
^ -0.0623 

When  these  values  were  run  through  RTRESP  as  initial 
conditions,  with  the  input  function  now  zero,  the 
results  were  inconsistent  with  common  sense.  For 
example,  the  response  of  the  system  with  an  initial 
condition  of  0.9938fps  of  forward  velocity  perturbation 
and  no  input  yielded  a change  in  U of  over  300fps  with- 
in ten  seconds  (the  initial  velocity  of  the  aircraft  was 
but  660fps ) . 

Using  the  state  transition  matrix  from  3ASMAT,  a 
program  called  ICRE3P  was  written  to  solve  the  response 
due  to  an  initial  condition  with  the  equation 

jc(t)  = (t)*x(0)  . 

The  results,  shown  in  Figure  6-7,  are  much  more  con- 
sistent with  the  earlier  figures.  In  the  output  por- 
tion of  the  print-out,  A1  refers  to  the  real  part  of 
the  first  root,  31  to  the  imaginary  part  of  the  first 
root,  and  likewise  for  A2  and  32.  The  matrix  E is  the 
matrix  of  coefficients  of  the  state  transition  matrix, 
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PL  PI 


PROGRAM  KRESP 


;i‘y 

30 
40 
50 
60 
?0 
SO 
90 
100 
110 
120 
130 
140 
ISO 
160 
170 

tso 

190 
2O0 
210 
220 
230 
240 
250 
260 
270 
2S0 
290 
300 
33  0 

w20 

330 

440 

350 

360 

370 


DIM  Et  16?  4 3?  DC  4?  1 3?XC4?4] 
MliT  READ  £ 


Dft T Hi  -5. 29. ’46-05?  O.  0035631975?  1.  1288176* -0.0015*091-73 
DATA  0.0100397* 1 ■ -0. 013779? -O. 0942425 
DATA  -0.0O01231O6? 1 . 87941E-06? 1 ? -0. 00069424 
DATft  2.  15602E-05?  8.001 13484* -0. 078526 * -5. 71775E-05 
DATA  -0.0OO212434? -O.00500833? 1 . 31718? -0. O0O442O4 3 
DATA  -0. 02OO49? 0. 13164849? 177.653?-0.  14041269 
DATA  -3 . 6054349E-05 ? -0 . 0058 1 48727 ? -0. 1 8 1 2 1 7 ? 8 . OOO 
DATA  -2. 09265E-05?  0. 080671698?  0. 224737 ? -0. 808219 
DATA  l?-8. 00856826? -1. 128878? 0. 08154636 
.DATA  -0.01 8039 ? - 7 . 9 1 1 35E-05 ? 0 . 0 1 062 1 36  ? 8 . 0942492 
DATA  8. 0O0123107? -1 . S706E-06? -8. 33737E-85?  0. 0086 
DATA  -2. 156E-05? -8.001 1849? 0.078517? 1 
DATA  -0 . 87 1 893 ? 0 . 68873556 ? -48 . 203644 ? -5 1 3 . 1 9324 
DATA  0.8O09086? -8.00605144? 0.39170456? 5. 1588733 
DATA  -7.49E-86?  7.51 33E-05? -0.004389544? -8.863124 
DATA  0. 0O196O4028? -8. 0001 15022? 0. 00427627*0. 0S29 
MAT  READ  D 

8. 9937396? -21 . 820388? -0. 035139? -O. 8623077 


9426 


95= 

68' 


TO 

TO 

TO 


4 

4 


DATA 

FOP  1=1  TO 
FOR  J=1  TO 
XC  I ? J 3=0 
NEXT  J 
NEXT  I 
FOR  1=1 
FOR  J=1 
FOR  K=1 
XC I ? J ]=EC 1+4* 
NEXT  K 
NEXT  J 
NEXT  I 
PRINT 
PRINT 
PRINT 
MAT  PRINT 
END 


J-4  ? K 3*DC  K?  1 3+XC  I ? J 3 


EXP  < A 1 *T  ) * 
C0S<B1*T) 


EXP<A1*T)  * 
SIN(B1*T> 


EXP<  A2*T> 
C0S<B2*T 


EXP 

31 


• A2  + T ■ 

n-:b2*t 


EXP< A1*T  > * 

cqs<bi*t> 


-8.226534030 

-21.80429874 

-0.035259895 

-0.023069852 


' EXP<  A1 *T)  * 
SIN<B1*T> 

13.062817347 

-10.21744129 

0. 133175189 

-O.022561 160 


EXP<A2*T)  * 
C0S(B2*T> 

1.220327219 

-0.814495920 

1 . 23026E-84 

-0.03923266O 


EXP<A2*r-  - 
SIN<B2*T  ■ 

20.03436953 

-0. 282251630 

2. 45808E-O3 

-3. 61228E-0-+ 


Figure  6-7 
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D is  the  matrix  of  initial  conditions,  and  X is  the 
matrix  of  coefficients  of  the  output.  The  time  his- 
tory plots  of  the  results  are  shown  in  Figure  6-8. 

Since  the  results  were  so  consistent  with  the  earlier 
program,  it  appeared  as  though  the  reason  these  figures 
were  not  coinciding  with  those  of  the  text  was  that  an 
error  existed  in  one  of  the  two  Melsa-Jones  programs. 


The  results  of  BAS MAT  (the  determinant,  the  in- 
verse, the  characteristic  polynomial,  the  eigenvalues, 
and  the  resolvent  matrix)  were  computed  by  hand  and 
found  to  be  correct.  The  prosram  RTRSSP  calculated  the 
same  eigenvalues  as  did  3A3MAT;  its  only  problem 
appeared  to  be  in  the  coefficient  of  the  sine  term  of 
the  pnugoid  root  in  the  forward  velocity  expression 
when  initial  conditions  are  present.  The  equations 
of  state  from  RTRE3P  were  verified  on  the  HP9830  by 
writing  a program  called  STRESP  which  calculated  the 
response  of  a system  to  a step  input  by  using  the  for- 
mula ( 

x(t)  = $(t)«x.(0)  + CaT  l*  ( t ) - 


where  [a]  is  the  plant  matrix,  ^b}  the  control  vector, 
k the  gain  (magnitude)  of  the  step  input,  and  $(t)  th 
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state  transition  matrix.  In  this  case,  x(0)  = 0. 

The  program  and  results  are  shown  in  Figure  6-9. 

Several  other  approximations  to  the  desired  input 
were  tried,  the  first  of  which  was 

S(,0  ~ T^T  ^ /Tt  J &*■  IT  «■  1 

Transforming  into  the  Laplace  domain, 

- f 

Jl(W]  ■ je 

The  solution  was  obtained  by  the  use  of  Bonnet's 


Theorem  (the  Second  Mean 

Value  theorem)  which  states 

that  rJ», 

i 

VS 
— » 

where  Gl(^)  s 

t 

f fLt)  It 

w 

For  -fit)  s e where  a. s and 

m --  f e' “<H 

= 1 - 

0 

-«tv 

and  jf'(t)  s - 

90 


fl  11 


PRGGRHM 


S i'RESP 


in  DIM  AC  4 * 4 ]*  FC  4 * 4 1 * QC  4*  4 1* RC  4*  4 J*  SC  4*  4 ] 

2n  DIM  BC  4*  4 ]*  CC  4 * 4 J * DC  4*  4 ]*  EC  4 * 4 ]*  FC  4*  1 ] 

3u  DIM  WC  4 * 1 3 * XC  4 * 1 3 * YC  4 * 1 3 * ZC  4 * 1 ] 

40  MAT  READ  B 

5n  DRTR  -5. 2974E-85* 0. 0085681975* 1 . 1288176* -8. 0015909672 
60  DRTR  0. 0100397 * 1 * -0.018779* -0.0942425 
70  DRTR  -0. 060123106* l . 87941E-06* 1 * -0. 00069424 
8H  DRTR  2.  15602E-05*  0.001 18484* -0.078526* -5.  7l77cr.E-05 
90  MRT  RERD  C 

1O0  DRTR  -0.0002 12434 *-0.00500833* 1.3 171 8 *-0.000442043 
10  DRTR  -0.020049*0. 13164849* 177. 653»-0. 14041269 
120  DRTR  -3.  6054349E-05* -0. 0058148727*  -0.  181217*  0. 00067.-::-: 
130  DRTR  -2. 09265E-05*  0.000671693*  0.224737* -0.  000219  Vr-: 
140  MRT  READ  D 

158  DRTR  1* -0.00856826.-1. 128373*0.00154686 

160  DRTR  -0. 010039*-?. 9tl35E-05*0. 01062186*0. 09424925 

1 70  DATA  0.000123107  * - 1 . 3706E-06* -8. 38737E-05*  9. 00069426 

180  DRTR  -2. 156E-05* -O.001 1849*0. 078317*1 

190  MRT  RERD  E 

200  DRTR  -0.071693* 0.60873556* -40. 203644* -513. 19824 
2.0  DRTR  0. 0009086* -0. 00605144* 0. 39170456*5. 1588783 
220  DATA  -7 . 49E-06, 7 . 5133E-05* -0. 804839544* -0. 063124955 
230  DRTR  0. 00 1960*02 3* -0. 000115022*0.00427627*0. 8829683 
240  MAT  RERD  R 

258  DRTR  8* -9. 64173* 638. 3179* 8136. 736 
268  DATA  0* -0.055395* -42. 6289 *-81.86226 
278  DRTR  8*0*8* 1 

288  DRTR  -0. 83 1056*8. 00290 l?»-0. 194487* -2. 45521 
290  MRT  READ  F 
300  DRTR  0* -69. 8* -26. 009* 0 
318  MAT  P*A*B 
.20  MRT  G*A*C 
:0  MRT  R=R*D 
340  MAT  S»A*E 
350  MAT  W*P*F 
360  MAT  W*<0. 02>*W 
370  MRT  X=Q*F 
380  MRT  X*<0. 82>*X 
390  MRT  Y»R*F 
4O0  MRT  Y*(0. 02>*Y 
410  MRT  Z=S*F 
420  MRT  2=<0.02>*Z 

430  PRINT  "THE  TIME  RESPONSE  OF  THE  STATE*  X(T)“ 

440  PRINT 

450  PRINT  "VECTOR  COEFF  OF  EXPC-2. 104333>T*CQS<3. 714?44>T 


FlCur*  6-9 
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460  PRIM  r UC  1 » l 3*  UC  2 * 1 It  ME  3»  1 ]» UC  4>  1 ] 

4 70  FRINT 

400  PRINT  "VECTOR  COEFF  OF  EXP(-2.  104338  1 T *SIN<3.  714744  • <" 
4'*0  PRINT  XCltl3tX[2tl3tXC3tl]tXC4tl3 
509  PRINT 

5.0  PRINT  "VECTOR  COEFF  OF  EXP  < -0 . 0045 1 1 029  > T*CGS  O . 062  7*  r 
5i.‘0  PRINT  VC  1 1 1 3»  YC 2 j 1 3>  YC  3t  1 3»  YC  4t  1 ] 

530  PRINT 

5<*0  PRINT  "VECTOR  COEFF  OF  EXPO0. 00451 1029>T*SIN<.0. 0627-47 
550  PRINT  ZC  1 * 1 3»  ZC  2 » 1 3*  ZC  3t  1 3»  ZC  4t  1 3 
560  END 


THE  TIME  RESPONSE  OF  THE  STATE » X<T) 

VECTOR  COEFF  OF  ZSF  2. 104338>T*COS<3. 714744>T 
0.233872332  \°  --312519  • 0.039193613  0.019483341 

VECTOR  COEFF  OF  t*P<  2. 104338)T*SIM<3. 714744)T 

-U. 045144590  1O.41547062  -0.117841333  0.021596520 

VECTOR  COEFF  OF  EXP . -0. 00451 1 029 >T*COS<0. 06274706 >T 
-318. 7879228  .206305933  -0.039138853  0.077583439 

VECTOR  COEFF  OF  EXP( -0. 00451 1029 >T*S IN <0. 06274706 > T 
-.•.3.35349032  0.075318325  -2.06336E-O3  -0.619O6O577 


« 


Figure  6-9  (continued) 
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Substitutin'?, 


Jilt)  d-t  = f(«-)  Gt")  - 


* ) At 


*rfO’e~tk)e  dt 


r^ci1  .it  2*  f '**  . 2jL  -<t  ->*  x. 

Jt  t cit  = 7 j e.  At  - * / e.  6 ^ 


lombinin.^  terms, 


.?«  A””-1" 


* -VJ? 

S 3 l* 


t+  --t 


_s PtF  j_ 

*4-5  i \ <*  *• 


J*  ■»  S 


For  T = 0.0 1 , <*  = “ 


3 Wife 


e e 


3 \H-U» 
tlSH  + * 


3 \vV  I U 
tl&Sl  -*  S 
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Multiplying  by  the  magnitude  of  the  input,  0.02, 


yielded 


0.02  <5  (t) 


6.2832 

62832  + 3 


The  magnitudes  of  the  results  of  this  input  were 
too  small,  the  probable  reason  being-  that  half  of  the 
function  exp^  /fv)  lies  to  the  left  of  the  t=0  axis 
and  is  contributing  nothing  to  the  incut. 

The  application  of  a 0.001  second  delay  to  this 
function  would  place  99^  of  it  to  the  right  of  the 
t=0  axis,  countering  the  problem  that  existed  originall 
A time  delay  of  0.001  seconds  corresponds  to  multiplyin 
the  Laplace  transform  by  0 ” 0 • 6 ;j , Approximating 
e“0.001s  by  the  series 

(0.00O  s (o.oo<)  S 

1 - o.oois  4*  — — ' Tv * • • ■ ) 
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R(s) 


becomes  j 

f0. OOOOOis2  - 0.0005s  + 0.05\ 

/ 6.2832  \ 

2 J 

\62832  + s ) 

= 79577.5 


- 500s  + 500,000 
s + 62832 


Ap-ain,  N(s)  is  of  greater  order  than  D(s). 


The  correct  answer  was  finally  reached  when  a 
variation  of  an  earlier  idea  (using  the  response  at 
t=l  from  a step  function  as  initial  conditions)  was 
pursued.  In  this  case,  the  two  steps  were  combined  into 
one  by  summing  two  step  functions,  one  commencing  at 
t=0,  the  other  at  t=l  and  negative  in  sign.  The 
e-raphical  representation  and  summation  is  shown  in 
Figure  6 — 14, 


1 

< 

P-**| 

i 

I 1 

1 

1 l 

-r  M 

'■•■T — 

i 

4 

-1 

1 

Figure  6-1  4 

In  the  time  domain,  the  result  is  written 

r(t)  = l(t)  - l(t  - t ),  where  r = 1. 
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The  Laplace  was  easily  found  and  was 

' i 

RCO  - i ~ s e 


Searching  for  an  approximation  of  e"s  which  would 
satisfy  the  requirement  for  the  order  of  the  numerator 
to  be  smaller  than  the  order  of  the  denominator  led  to1 


\-\+ 

J + 1*  £«•  ^+  - 

4S  - an*  -*•  <o  - *7  •*!: 
ws  + ti  + s 


Substituting  into  R(s)  = 0.02(--je."  ) yielded 

x „ r s'- 

= o.ohi— 7-7 ? 

-*•  * Ls  + s* 


The  results  of  the  RTRESP  computations  and  the  time 
history  plots  are  shown  in  Figures  6-15  and  6-1 6.  Their 
agreement  with  the  text  figures  confirmed  the  accuracy 
of  the  RTRESP  (for  other  than  initial  condition  inputs) 
and  GRAPH  programs  and  provided  a valid  approximation 
of  a finite  impulse-type  input. 


''•-Ogata,  Modern  Control  Engineering.  Prentice-Hall, 
1970, 
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VII.  LONGITUDINAL  FEEDBACK 


The  primary  purpose  in  usingr  feedback  in  a control 
system  is  to  reduce  the  system’s  sensitivity  to  external 
perturbations.  In  the  open  loop  system,  the  area  of 
main  concern  is  calibration?  in  the  closed  loop  sys- 
tem, it  is  stability.  Use  of  feedback  can  enhance  the 
rapidity  of  response,  reduce  errors  of  response,  and  in- 
crease damping  to  satisfactory  levels;  it  can  also  over- 
correct errors  or  decrease  damping  to  the  point  where  an 
otherwise  stable  system  becomes  unstable.  It  is  of  in- 
terest, then,  to  be  able  to  "measure"  the  stability  of 
a system. 

The  stability  of  a system  may  be  ascertained  by 
an  analysis  of  the  characteristic  equation  of  the  system, 
that  is,  the  denominator  of  its  transfer  function. 

Simply  stated,  if  the  real  part  of  every  root  of  the 
equation  is  negative  in  sign,  the  system  is  stable. 

That  this  is  so  is  manifested  by  the  relationship  between 
the  roots  of  the  characteristic  equation  and  the  solution 
of  the  differential  equation  which  it  represents. 
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For  a multiple  input/output  system,  the  feedback 
effects  on  the  transfer  functions  must  be  dealt  with  by 
matrix  manipulation 

The  solution  to  the  system  of  differential  equations 
x = \_A]x  + {b}u 

was  see"  to  be,  in  the  Laplace  domain, 

x(s)  = [sI-a]  (b]u(s)  = [<$>  (3)]{b}u(s). 

The  ope"  loop  transfer  function  was  then 

-»'■>  * 5$  - 

where  cj>(s)  is  the  resolvent  matrix.  In  a 4x4  system, 
the  transfer  function  is  written 
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% 

Z ^ /a 

tv 

l - y 

Z Vi  /a 

Gq ( s ) - 

M5!.  <P„  ft* 

^ 'fnv  «♦»■»  <fW 

, J 

- \ 

2<<ViVi/  A 
L 2 

A 


j=  4 


The  last  quantity  will  be  written 


where  A = d e t I - a] 


Applying  the  feedback  signal  to  a particular  parameter, 
for  example,  Xj,  yields 


H = [0  0 h O] 


where  the  block  diagram  is 


(The  solid  arrows  are  used  to  signify  a multiple  input/ 
output  system.)  The  closed  loop  transfer  function  is 
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n(  g ) - ^Q^S) 

J'S)  - 1 + GnK 


In  matrix  format,  ^G(s)J  = ^1  + G0h1  |g0(s)|. 

To  find  this  matrix  product,  the  matrices  were  taken 
term  by  terms  r 

1 0 hGi  0 

rT  „ 0 1 hG2  0 

j_I  + G0Hj  - o 0 l+hG3  0 

0 0 hGi;  1 


[i  + g0h]“  = 


1+hC-o 

0 

0 

0 


0 -hGi  0 

l+hG^  -hG2  0 
0 1 0 
0 —hG/;  1 +hG^ 


/here  = det  [i  + G0h]  = 1 + hG^. 


6, 

&». 

~cT  ~ 


T bj  'P.j 

2 b> 

2 bj  'fxj 

. 2 bj  ^ 


2 bj  ^ 

21  ^ 

2 Vj 


(\  * W 


2T  b;  ^ ' 

j 2 bj  <f»j 

' 

2 bj  'f’tj  , 
A+  K"Z  bj  4>3j 
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The  closed  loop  transfer  function  is  ther 


( -v  -v  «f,,  -+  ^ 

) V>, vfv«  -*  ^ -t  \*3  ^)ki  ^ 

"*  -v  r ’'Pit 

^ ^1  'fw  ■*  w V,  v£,,  y-  \*V4 

& + V\(>,ifci  + W'fn  t-  * 'omf,*) 


In  the  general  form,  the  first  subscript  of  »4>  in  the 
denominator  will  vary  as  the  feedback  parameter  is 
varied.  Regardless  of  the  form  or  magnitude  of  h,  only 
the  denominator  is  affected;  the  numerator  is  a function 
solely  of  the  openloop  resolvent  matrix,  $(s),  and  the 
control  vector,  {^b^. 

The  most  important  characteristic  of  the  transfer 
function  is  that  the  frequency  response  may  be  obtained 
by  letting  s = ju>  , where  is  the  applied  sinusoidal 
forcing  function.  That  this  is  so  is  shown  by  assuming 
that  the  output  contains  a sinusoidal  component  of  the 
same  frequency  as  the  input  in  the  form 

A*  R(ti  ) • sin  t + (u  ) j 


where  A is  a constant,  R(o)  is  the  magnitude,  and  c$(ui) 
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is  the  phase  angle  between  the  input  and  the  response. 
For  an  input  of  Asini^t,  the  Laplace  of  which  is 

Acs 

S2  + ta,1- 


the  output  is 


= G(s)r(s)  = 


At4  G ( s ) _ A co  G ( s ) 

s^  + wj"  (s  + j«>)(s  - jui) 


Using  partial  fraction  expansion: 


y(s)  ’ i ...a 


+ transient  terms 


Assuming  the  transient  solution  is  stable  allows  the 
transient  terms  to  be  dropped,  so 

n Aul  G-U)  - 

Jjx  - ' ” r ■ — 


Aui  6CO 


i + * _ 

J l - ♦O'J 


A & (jui) 


Solving  for  ,y(t)  by  taking  the  inverse  Laplace  transform: 


y(t)  = 


-.at  A6*(nv*)) 

— - — 4.  -I /-v  * 


e + 


(7-1) 


Since  G(jiA  ) is  complex,  write  it  as 


;(jU)  = x(yi)  + jY(ui)  and  G(-joi)  = x(ui)  - jY(v^  ) 
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where  | G(  jo)  | = + Y2(o)  and  argG(ji»  = tan_1^| 

Substitute  into  equation  (7-1): 

3W  = ' , j Y(j)]e*J'iC 


AXW^^1)  + e^cJ 

>1 


= AXU)5'*^  + AYU)  c*>  ut 


where  x = Yaw 


->  vw 

1 XU) 


This  is  the  same  form  as  the  assumed  input,  so  that 
R(s)  becomes  J X2(ui)  + y2(u})  = | G ( j ui)|  and  the  phase 

an£rle<f(oi)  becomes  tan"  - Y(u>  )/X(^)  = areG(  jui) . 

It  becomes  a simple  matter  to  apply  a feedback  to 
a plant  matrix,  obtain  the  transfer  functions  in  each  of 
the  variables,  and  plot  the  magnitude  of  the  responses 
to  a sinusoidal  input  against  the  frequency  of  that 
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input 


o This  plot  is  termed  a Bode  (Eo’-dee)  plot,  named 
for  its  developer,  H.W.  Bode,  One  of  the  important 
uses  of  the  plot  is  in  measuring-  the  stability  of  e 
system. 


In  order  to  facilitate  construction  of  a Eode  plot, 
and  in  order  to  avoid  relying  upon  asymptotic  approxi- 
mations of  magnitude  and  phase  angle,  a program  named 
BODE  was  written  for  the  HP9330  computer/plotter.  Using, 
as  input,  the  stability  derivatives  and  the  type  and  mag- 
nitude of  feedback,  the  program  prints  out  the  resulting 
transfer  functions  (in  the  Laplace  domain),  repeats  the 
feedback  inputs,  and,  for  given  limits  of  frea_uency  and 
magnitude,  constructs  the  3ode  diagram  on  semi-log  paper. 
The  feedback  form  may  be  linear,  integration  (i/s),  or 
diff erentiation  (s). 

The  program,  listed  in  Figure  7-2,  is  explained: 

line  10  dimensions  open  and  closed  loop  matrices 

11  20-60  defines  variables  in  terms  of  stability 
derivatives ; 

XrJ  = Xu 
X’Y  = Xw 
XDE  = Xjt 

MWDOT  = , etc. 

line  70  iterates  program  for  x-i  , X2.  x^,  and  X4 
(u,w,q, 8)  outputs 
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H i'.  i i CM 


Ji  i *l«t 


22iiov7© 


I Vj 
hi  *J 
30 

4 8 
50 
6y 
70 

90 
100 
110 
120 
130 
14t;1 
150 
160 
170 
1 30 
190 
200 
1 7 

. c.O 

- r0 


ilfb.O  I,  :,L 

. i ’ ’ r.  J 

.71  , Z2  * 23 
M 1 * M2 ? M3 
M4 . M5  > M6 
>.4  IS  X Q 

J=1  TO  4 
15  , 30.- 
45"*" 


III' 2 XO.XW.XDE 
liKE  ZO,  ZW,  ZDE 
HRE  hu> nw> MM HOT 
OPE  MR, MQ?  MDE 


290 
SOU 
3 x 0 
320 
330 
240 
350 
3 S 0 
370 
3dU 
390 
4O0 
410 
420 
430 
440 
450 
460 
470 
480 
490 
500 


IHH  HI  8 , o 8 1 

KLl'i 
PE  H 
REM 
REM 
REM 
FOR: 

WR I IE 
FORMAT 
PRINT 
DEG 
T3=y 

00-1016.9 
G- 32. 073 
X 1 =-0 . 0097 
X2=0. 0016 
X 3=0 
X4=0 

21--0. 0955 
’2=-l . 43 
"3 --69.3 
; l ! =- Ti 
17-  0235 

r,l  00  1 3 

M4 =00*112 
i'!5=- 1 . 92 
1 

Or-  1 = 1 TO  7 

uf  n-o 

FOR  r-.  = 1 TO  7 
1IC  1.10=0 
NEXT  K 
NEXT  I 
DC  5 1=1 

NC  3. 4 3=M6+23*M3 
NC  4,2  1=23 
NC  4,  1 ]=X3 
DL4  1=-M5-Z2-X1 

HI  2,4  3 = X 3 * •:. Z 1 * M 3 + M 1 > +Z  3 * < M 2 - X 1 * M 3 •' - M6 * ( X 1 f22> 

NC 3 , 2 3=X3*Z 1 -23* < X 1 +M5 > +M6*U0 
HC 3, 1 J=-X3* ( Z2+M5  > +Z3*X2+X4*M6 
DC  3 ]=Z2*M5-X2*2 1 *X 1 * < M5+Z2 > -M4-M 1 *X4 
HC 1 ,4 ]=X3*<Z1 *M2-Z2*M1 )+Z3*<Mi*X2-M2*Xl >+M6*<Z2*Xl 
NC  2 , 2 J=X3* < US*M 1 -Z 1 *M5  > +Z3*X 1 *H5-U0*M6*X 1 +X4* < M6*Z 1 -H i *Z3 > 

NC  £ , 1 ]=X-3* l-  Z2*M5-M4.)  -Z3*<X2*M5+G*M3  >-M6*G+M6*X2*U0-M6*Z2*X4 
DC  2 3=-X  1 * Z2*M5-M4.> +Z 1 * < X2*M5+G*M3 > -M 1 * i X2*U6~G  > +X4*  < M 1 *Z2  - ,12 *Z  i > 
NC  5, 4 3=NC  5, 3 3=0 
HC  4,3  ]=M6 

NC  3 , 3 3=-M6*X 1 -M6 *Z2+M 1 *X3+M2*Z3 

NC  2,3 3=M6*X 1 *Z2+M2*X3*Z 1 +M1 *X2*Z3-M 1 *X3*Z2-M6*X2*Z1 -M8*X 1 *23 


Figure  7-2 
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■ • ! »!  i '[  l i 3 1=0 
■ I ! HI  l • 2 1=G:»  < -til  -i’lt>*z  i 
5 !m  HC  1*1  ]=G*  '.116  O:  Z3*M2> 

540  o::  i 3 - g • z l * n3  - ti i + Z2 > 

550  DISP  IS  FDBK  0*  LINEAR*  5*  OR  1V‘5 

5 GO  INPUT  Z 

5,-0  IF  Z=0  THEN  940 

580  IF  Z=1  THEN  610 

590  IF  2*2  THEN  680 

600  IF  2=3  THEN  790 

610  DISP  "INPUT  FEEDBACK  GAIN"? 

620  1 1 IPUT  K 

630  PRINT  "FEEDBACK  IS  LINEAR" 

640  FOR  1=1  TO  7 
650  DL  I J=DC  I 3-NC  I > J ]*K 
660  NEXT  I 
670  GOTO  960 

680  DISP  "INPUT  FEEDBACK  GAIN"? 

F90  INPUT  K 

" PRINT  "FEEDBACK  cORM  IS  K*S" 

10  for  1=7  TO  1 slip  -1 
"So  SC  I + 1 * J 3=NC  I ? J 3*K 
" 30  NEXT  I 
740  SC  1 > J 3=0 
"50  FOR  1=1  TO  7 
7 60  DC  I 3=DC  I 3-’3C  I * J 3 
770  NEXT  I 
730  GOTO  960 

79o  DISP  "INPUT  FEEBACK  GAIN"? 

800  INPUT  K 

810  PRINT  "FEEDBACK  FORM  IS  K/S" 

820  FOR  1=7  TO  1 STEP  -1 

839  DC  1 + 1 ]=DC  I 3 

840  NEXT  I 
850  DC  1 3=0 

86©  FOR  1=1  TO  7 
870  DC  I ]=DC  I 3-NL  I * J J*k 
880  NEXT  I 

890  FOR  1=7  TO  1 STEP  -1 
900  NC 1 + 1 * J 3=NC I > J ] 

9lo  NEXT  I 
920  NC  1 » J 3=0 
930  GOTO  960 

940  PRINT  "FEEDBACK  IS  NOT  APPLIED" 

950  GOTO  970 

960  PRINT  "FEEDBACK  GAIN  IS  "K 

970  PRINT  "TRANSFER  FUNCTION  COEFFICIENTS  C INCREASING  P ONERS  OF 
980  PRINT 

990  PRINT  "COEFFICIENTS  OF  NUMERATOR" 

1000  FOR  1=1  TO  7 


Figure  7*2  (continued) 
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U1 L Li  Kim  Ml  1 J J 3 

1020  next  l 

l 'J  >y  FRIMT 

1040  PKlNi  COEFFICIENTS  OF  DENGMINhTuR  • UlhR  FJ.4 
10 sy  FOR  1=1  TO  7 
1000  FRIMT  DII3 
1070  ME XT  1 

1080  LISP  "PLACE  PAPER  OH  PLOTTER  5 
1O90  STOP 

11O0  LISP  " INPUT  G/  MAX  < LB > " 5 
1110  INPUT  G9 

1120  Ii ISP  "INPUT  OMEGA  MINIMUM"; 

113©  INPUT  00 
1140  GOSUB  1650 

1150  FOR  L*LGT<O0)+1  TO  LGT (OB) +5 
1160  .:=10TL 

1170  FOR  W=X/10  TO  X STEP  X/50 

1 1 30  XT  1 3=  ( NT  1 5 J 3— NT  3 * J 3*Wt2+NT  5 , J ]*W  i 4 -NT  7 , J 3* U f 
1 1 90  XL  J 3= XT  J ]+< NT  2 > J 3+W-NC  4 » J ]*Wt3+Nl  6 » J J*W*5  > *Z 
12O0  XT  J 3=SQR<XC J 3 > 

I2l0  VI  J ]=■:  DT  1 3— DT  3 J*Wt2+Bt  5 3*Wt4>  t2 
1 220  YT  J 3=YC  J 3+  < DT  2 3*I4-DC  4 3*W  1 3+BT  6 3 Hi  r 5 1 *2 
1230  YT  J 3=SQR< YT J 3 > 

1 240  RT  J 3= XT  J 3-  YT  I 3 

1750  IF  20*LGT- RT  J3KG9-140  THEN  1270 
1 260  PLOT  LGT <w>  * 20*LGT  < RT  J 3 ) 

1270  NEXT  W 
1280  NEXT  L 
1290  PEN 

1 500  PLO T LGT  < 00  ) +0 . I > G9- 1 > 1 

1310  LABEL  1.2> 1.7»0»7/10> “DECIBELS" 

1320  PLOT  LGT < 00 > +2. 2>  G9-140J 1 
1 330  LABEL  < * > 1 . 2 ->  1 . 7 > 0 * 7/ 1 0 > “ RAD  I ANS  'SECOND  " 

1340  PEN 
1350  GOSUB  1910 

1360  FOR  L=LGT'.O0>  + 1 TO  LGT<O0>+5 
1370  X*10tL 

1580  FOR  W-X/10  TO  X STEP  X/50 
1390  T9=T3 

1400  Q 1 =NC  2 » J 3*W-NC  4»  J 3*Wt3+N[ 6*  J 3*14+5 

1410  Q2»NC 1 j J 3 -NT  3? J 3*W+2+NC  5>  J 3*W+4-NT  7<  J 3*H r 6 

1 420  Q3*DT 2 3*W-DT  4 3*W+3+DT 6 3*W+5 

1430  Q4=DC  1 3-BT3  3*W*2+DC5  3*Wt4-DC7  3*WT6 

1440  GOSUB  2030 

1450  T3=T 1 -T2 

1460  IF  ABSCT3-T9K180  THEN  1500 
1470  T3=T3+360 

1480  IF  ABS<T3-T9X1S0  THEN  1500 

1490  T3*T3-720 

1500  IF  T3C-300  THEN  1520 


Figure  7-2  (continued) 
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I 5 ' u PLOT  L G I * i'4  1 ’ 7 J 
t -.LO  NEXT  u 
l.'U  NEXT  I 
1040  PEN 

1350  DISP  IHPlJl  LETTER  SIZE  < 

1560  INPUT  Z 

1 5 1 0 LABEL  \ + ? Z?  2 ? O ? 7 1 0 > 

1 5 SO  DISP  “YOU  ARE  IN  THE  LETTER  MODE 
1590  LETTER 

1 DO  DISP  “ANY  OTHER  LETTERING" 5 

1610  INPUT  Z 

1620  IF  Z=i  THEN  1550 

1630  NEXT 

1640  STOP 

1 658  SCALE  LGT  ( GO  > -0 . 5 ? LGT < 00  .'+5.5’  G 9 - 145?  05  + 
1660  DISP  "CHECK  ALIGNMENT  OF  PLOTTER" 

1670  PLOT  LGT(O0)?G9-140? 1 
1680  STOP 

1 690  PLO T LOT ( 00  > +5  ? G9 ? 1 

1700  DISP  "RF -CHECK  ALIGNMENT?" 5 

1710  INPUT  1 

1 720  IF  Z=0  HEN  1740 

1730  GOTO  1660 

1740  VAX  1 3 LGT  < 00 >?20» G9- 140?G9 
1750  LABEL  ( + ? 1? 1.7»0,7/10> 

1760  FOR  Y=G?-140  TO  G9  STEP  SO 

1770  PLOT  LGT(O0)?Y?1 

1 780  CP LOT  -5? -0. 3 

t',90  LABEL  ■ 1 800 ) Y 

1 8O0  F ORMAT  F4 . 0 

1810  NEXT  V 

1820  FOR  X=Lj ! <00  • + 1 TO  LGT (00 > +5 

1 830  PLOT  X « G9- 140? 1 

1340  CPLOT  -4. 5, -1.5 

1350  LABEL  1860  -1 0tX 

I860  FORMAT  F3.2 

1870  NEXT 

1830  PEN 

1890  RETURN 

1900  STOP 

1910  SCALE  LGT  ( 00  > -0.5?  LGT  < 00  > +5.5? -325  ? 425 
1920  DISP  “RE-CHECK  ALIGNMENT  OF  PLOTTER"? 
1930  PLOT  LGT (00) ? -300? 1 
1940  STOP 

1950  PLOT  LGT (00) +5? 400? 1 
I960  STOP 

1970  LABEL  ( *?  1 . 2?  1 . 7?  O?  7-  10) 

1980  FOR  Y=-200  TO  100  STEP  100 
1990  PLOT  LGT <00 • +4. 7? Y» 1 
2000  CPLOT  O- -0.3 


Figure  7*2  (continued) 
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-'  .i  i h L H b F L <:  2050  , 

20  20  ML' XT  V 

2 0 3 0 PLOT  L G T u 0 > f A . 9 5 > 9 5 , 1 

2040  LABEL  <*>  1 . 2 > l . 7 j 0 > 7 ■ 1 0 > " DEG 

2 0 5 0 F 0 P M Ft  T F 4 . lit 

2060  RETURN 

2078  STOP 

20 SO  IF  01  0 HUD  02 <0  THEN  2130 

2000  IF  01 <0  HMD  02< 0 THEN  2150 

2100  IF  Ol<0  AMD  02 >0  THEN  2170 

2110  T 1 = ATN  < 0 1 ’02 >-180 
2120  GOTO  2180 
2130  1 l=FtTIKQl/Q2> 

2140  GOTO  2180 
2150  T 1 =ATN < 0 1 '02 > 

2160  GOTO  2180 
7 1 70  T 1 - ' s0+p  ! N i,  Q 1 /Q2  / 

7 1 30  If-  ■ :• ; 0 HMD  04<O  THEN  2230 

. ! IF  03x0  RND  04 <0  THEN  2250 

.00  IF  03< 0 RND  04 >0  THEN  2270 

..  10  T2=ATN<Q3/Q4> 

-220  RETURN 

2230  T2*18Q+AT N < 0 3 / Q 4 

2240  RETURN 

2250  T 2=  1 30+ATN*. 03-  04  > 

2260  RETURN 

2270  T 2=360+RTN  < 03/04  ) 

2280  RETURN  / 

2230  END  s 


Figure  7-2  (continued) 
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11  80-110  format 


11  1 20-270  assigns  values  to  variables,  initializes 
T3  (phase  angle,  used  for  comparison) 

11  280-330  initializes  numerator  and  denominator 
coefficients  to  zero 

11  3^0-540  computes  values  of  numerator  and  denomi- 
nator coefficients 

11  550-1070  inputs  feedback  type  and  magnitude,  com- 
putes and  prints  closed  loop  numerator  and 
denominator  coefficients 

11  1080-1640  inputs  frequency  and  magnitude  limits 
(of  response),  labels  graph,  and  plots 
magnitude  and  phase  angle  against  applied 
frequency 

11  1 650-1900  subroutine  for  scaling  and  labelling 
magnitude  graph 

11  ’p1 0-2070  subroutine  for  scaling  and  labelling 
phase  angle  graph 

11  2080-2280  subroutine  which  computes  phase  angle 
correct  for  sign  and  quadrant 

line  2290  EA'D  statement 


To  verify  the  program,  the  McRuer  example  was  again 
used.  The  text  results  are  shown  in  Figure  7-3»  The 
BODE  results  are  shown  in  Figures  7-4,  7-5,  and  7-6  for 
the  three  parameters  common  to  both  programs,  * , and 


The  comparison  of  the  ot  values  of  the  text  and  the 
graph  of  w from  BODE  needs  some  explanation.  The  text's 
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magnitude  for  of  at  low  frequency  is  ~4.0;  the  results 
of  BODS  showed  the  magnitude  of  w at  low  frequency  to  be 
-*  61  , a 57  decibel  difference.  For  U0  = 660, 

w = otU0  = 660« 

log(w)  = log(660)  + log  («.). 

Since  log(660)  = 56.4,  the  57  db  difference  is  correct.  • 

Table  7-2  compares  the  results  of  the  BODE  program's 
transfer  function  coefficients  with  those  of  the  text. 

NUMERATORS 

BODS  0.?u68sI 2 *  - 791 .5l43s  - 1144.4512 
u 

text  0.1103s2  - 799s  - 1151 

BODS  69.8s3  + 26675.783s2  + 258.7485s  + 79.94357 

W 

text  69. 8s3  + 17350s2  + 168.6s  + 80.2 

BODE  26.00926s2  + 35.93499s  + 0.35011 

6 , ? . 

text  26.12s2  + 36.1s  + 0.351 

I DENOMINATOR 



BODE  s4  + 3.3597s3  + 26.6754s2  + 0.26271s  + 0.07198 
text  s4  + 4.22s3  + 18.28s2  + 0.2097s  + 0.0724 

Table  7-2 
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Inasmuch  as  the  text  figures  aopeared  to  be  slide-rule 
accuracy,  the  results  of  EODE  were  considered  more 
accurate.  The  BODE  graphs,  which  were  computed  with 
zero  feedback,  compared  very  favorably  with  the  text, 
and  the  program  may  be  considered  verified. 


The  transfer  functions  may  also  be  written  as 
functions  of  the  stability  derivatives,  and,  in  this 
way,  the  effects  of  various  types  of  feedback  may  be 
predicted.  Writing  the  equations  of  motion  (equations 
5-1 ) in  the  Laplace  domains 
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where  A 
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The  transfer  functions  for  w(s)/J<»(s),  q(s)/Jc(s),  and 
£(s)/J«  (s)  are  similarly  found.  Compiling-  the  coef- 
ficients for  each  of  the  numerators  yields  Table  7-3* 


The  time  history  response  may  be  obtained  by  sub- 
stituting the  error  signal  for  the  input  term  to  the 
riant  matrix.  In  other  words,  where,  for  the  open  looo 


svs 'em , 


£(t)  = + ^b}u, 


in  the  closed  loop  system 


*(t) 


x(  t) 


u = r - kx,  and  x(t)  = [Ajx  + {b}  (r  - |.k]{xj)  . 


Rewriting, 


x(t)  = (A  - bk^)x  + br. 


UMERATORS 
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As  originally  formulated, 


W“u  "w 
i»  u Mw 


Zw  'o+Zq 


H ■ 


Assuming  the  input  to  be  zero 
(r  = 0),  the  effects  of  a spe? 
error  feedback  can  be  seen  by 
letting  |_kj  = l*u  0 0 oj 


and  solving  for  - bk * ^ : 


'fat.  fen 


o o 


Xu  - X<e,Wu 

ZM  - L<rt  k„ 

Mu  - h* 
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This  is  the  new  plant  matrix  from  which  the  eigenvalues 
may  be  calculated  and  roots  plotted.  As  ku  is  varied,  the 
roots  move,  and  this  movement  is  the  critical  factor  in 
retaining,  gaining,  or  improving  stability.  Plotting 
the  roots  also  shows  how  variations  in  k affect  the 
undamped  natural  frequencies  and  damping  radios  and  gives 
a rough  idea  of  the  sensitivity  of  the  system. 

In  order  that  the  roots  of  a given  system  subject 
to  feedback  could  be  quickly  ascertained,  a program 
named  FDFK  was  written  for  the  HP9830.  Using  as  inputs 
the  resolvent  matrix  of  the  open  loop  system,  the  control 
Vector  {b},  the  coefficients  of  the  open  loop  character- 
istic equation,  and  the  type  and  magnitude  (gain)  of  the 
feedback,  the  program  computes  [_A  - bkj]  , finds  and  prints 
the  new  roots,  and  calculates  the  resulting  frequencies, 
periods,  and  damping  ratios.  See  Figure  7-7  for  a listing 
of  FDBK.  The  REMark  statements  in  the  listing  provide  an 
explanation  of  the  procedures  used. 

Figure  7-8  is  an  example  of  the  program's  output. 

In  this  case,  the  McRuer  example  was  again  used  with  no 
feedback  applied;  the  results  were  in  agreement  with  the 
examole . 
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PROGRAM  FDBK  08DEC76 


10  REM  THIS  PROGRAM  TAKES  A 4X4  RESOLVENT  MATRIX  FROM  DATA 
20  REM  STATEMENTS?  CONSTRUCTS  THE  OPEN-LOOP  TRANSFER  F'lMi:rHTN? 
30  REM  AND  FINALLY,  BY  APPLYING  FEEDBACK  FROM  ANY  ONE  OF  f HE 
40  REM  FOUR  PARAMETERS,  CONSTRUCTS  THE  CLOSED-LOOP  TRANSFER 
50  REM  MATRIX?  THE  FEEDBACK  MAY  BE  LINEAR,  K-S,  K*S*  OR  MOT 
60  REM  APPLIED  AT  ALL.  ANY  VALUE  OF  FEEDBACK  GAIN  MAY  PE  USED. 
70  DIM  At4,4I,BC4?4],CC4,41,DC4,43,EC4,4],FC4,4]?MC4] 

80  REM  INITIALIZE  COEFFS  TO  ZERO 
96  FOR  1=1  TO  4 
10O  FOR  J=i  TO  4 

110  AC  I,  Jl-BC  I,  J3*CE  I,  J3*Ot  I,  J3-EC  I,  J]=FC  I , J 3=0 
120  NEXT  J 

130  GC  I 3=H[  I 3=  It  I 3=J[  I 3*KC  I 3-HE  I 3=0 
140  NEXT  I 

1 50  S = T = U = V = W = X = Y = X 9 = Y 9 * 0 
160  PRINT 

170  WRITE  <15,180) 

1S0  FORMAT  45"*" 

190  FORMAT  22"  *" 

2O0  PRINT 

210  REM  ENTER  RESOLVENT  MATRIX  IN  INCREASING  POWERS  OF  S 
220  MAT  READ  B 

230  DATA  0,0.7567,-45.046,-627.429 
240  DATA  0,0.004061, 3.0646,5.91353 
2 50  DATA  0,0.0. -0.07797 
-60  DATA  0.00242,-0.00023,0.01401,0.1893 
270  MAT  READ  C 

280  DATA  19.4854,0.004443,-31.144,-135.5619 
290  DATA  -0. IS 365, 0. 02695, 6. 401356, 3. 075099 
300  DATA  0.00242,-0.00023,0.01401,-0.0039928 
310  DATA  0.000124,-0.0235,1.4397,19.52637 
320  MAT  READ  D 

330  DATA  4.21  0.0016.0,-32.2 
240  DATA  -0.0355,2.  7897,660,0 
350  DATA  0.000124. -2.0235, 1.4397,0 
360  DATA  0,0»  * >4,2;  ?? 

370  MAT  READ  E 
330  DATA  1,0, 0,0 
390  DATA  0, 1,0,0 
400  DATA  0,0, 1,0 
410  DATA  0,0,8, 1 

420  REM  ENTER  CONTROL  VECTOR  <B> 

430  MAT  READ  M 

440  DATA  0,-69.3,-26.009,0 

450  REM  ENTER  COEFFICIENTS  OF  ORIGINAL  CHARACTERISTIC  EQUATION 
460  REM  IN  INCREASING  POWERS  OF  S 
470  READ  N0,N1,N2,N3,N4 

480  DATA  0.07797-0. 19329, 19. 52637, 4. 2197, 1 
490  DISP  "IS  FDBK  0-  LINEAR,  K*S,  OR  K^S"? 

500  INPUT  L 


Figure  7-7 
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510  DIsP  -WHICH  PARAMETER  DESIRED  FOR  FDBK  " • 

520  INPUT  Z 

5 30  DISP  "ENTER  OR  IN " 5 

540  INPUT  K 

550  IF  Z#1  THEN  570 

560  PRINT  FEEDBACK  PARAMATER  IS  HORIZ  VELOCITY” 

570  IF  Z*2  THEN  590 

580  PRINT  'FEEDBRCK  PARAMETER  IS  VERT  VELOCITY" 

590  IF  Zi3  THEN  610 

6O0  PRINT  "FEEDBRCK  PARAMETER  IS  PITCH  RATE" 

610  IF  Z#4  THEN  630 

6Z0  PRINT  "FEEDBACK  PARAMETER  IS  PITCH  ANGLE” 

630  PRINT  "GAIN  = "K 

640  REM  MULTIPLY  PHI <S>  TIMES  B (.CONTROL  VECTOR  - TO  FORM 
650  REM  NUMERATOR  OF  GO < S > (OPEN  LOOP  TX  FN> 

660  FOR  1=1  TO  4 
670  FOR  J=1  TO  4 
6S0  HC  I 3=B£  I . J 3*MC  J ]+Ht  I ] 

690  IC  I ]=CC  I*  J]*MCJ3+IC  13 
'O-  1C  I 3=DC  1 1 J ]*MC  J ]+JC  I I 
’10  KC  I 3=EC  I « J 3*M[  J ]+KC  I ] 

720  NEXT  J 
730  NEXT  I 

740  REM  DENOMINATOR  OF  G0(S>  IS  THE  CHARACTERISTIC  EON  -.SEE 
750  REM  LINES  440-470).  NUMERATOR  OF  CLOSED-LOOP  TRANSFER 
,'60  REM  MATRIX  IS  THE  SAME  AS  THAT  OF  THE  OPEN-LOOP  MATRIX. 
770  REM  THE  DENOMINATOR  IS  OF  THE  FORM  H TIMES  SUM<BCJ)*PHI 
780  REM  <Z« J PLUS  CHAR  EQN  OF  PLANT  MATRIX. 

790  C0P  J=1  Ti;  4 

8O0  1 =BL  Z>  J 3*M[ J 3+T 

310  U=CCZ«  J3*M[  J3+U 

820  V=DC  Z f J 3*MC J 3+V 

330  W=ECZ>  J!*MC  J3HI 

340  X*F[  Z • .*  3*MC  J 3+X 

850  NEXT  J 

360  IF  L=0  THEN  990 

870  IF  L=1  THEN  1030 

880  IF  L=2  THEN  1110 

390  IF  L=3  then  9O0 

90O  PRINT  "FEEDBACK  MODE  IS  K/S” 

910  PRINT 
920  S=T*K 
930  T»U*K 
940  U=V*K. 

950  V*W*K 
960  W=X*K 
970  X»0 
980  GOTO  1200 

990  PRINT  "FEEDBACK  IS  NOT  APPLIED" 

1000  PRINT 
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luio  t-u=v=w=x-y=o 

1020  GOTO  1200 

1030  PRINT  "FEEDBACK  MODE  IS  LINEAR" 

104O  PRINT 
1050  T=T*K 
1060  U=U*K 
1070  V=V*K 
1 030  W=W*K 
1090  X=X*K 
1100  GOTO  1200 

1110  PRINT  "FEEDBACK  MODE  IS  K*S" 

1120  PRINT 
1130  Y=X*K 
1140  X=W*K 
1150  W=V*K 
1160  V=U*K 
1170  U=T*K 
1180  T*0 
1190  GOTO  120O 

1200  REM  SUMMING  TERMS  IN  DENOMINATOR  OF  CLOSED-LOOP  TX  MATRIX 

1210  S=S 

1220  T=T+N0 

1230  U=U+N1 

1240  V=V+N2 

1250  W=W+N3 

1260  X=X+N4 

1270  Y=Y 

1280  WRITE  <15* 1 90 > 

1290  PRINT 

1300  PRINT  "THE  COEFFICIENTS  OF  THE  NUMERATOR  OF  THE  CLOSED-LOOP  r 
1310  PRINT  " - FROM  S*<0>  TO  ST<3> > " 

1320  PRINT 

1330  FOR  1=1  TO  4 

1340  PRINT  "COEFFS  UF  G<"I")M 

1350  WRITE  <15- *470>HtI]»  IC  I ]>  JC  I 3»KC  II 

1360  PRINT 

1370  NEXT  I 

1380  WRITE  <15*190) 

1390  PRINT 

1400  PRINT  "COEFFICIENTS  OF  THE  DENOMINATOR  ARE:" 

1410  PRINT  " (FROM  St<-1>  TO  St<5>>" 

1420  PRINT 

1430  WRITE  <15* 1470>S*T,U,V 
1440  WRITE  <15* 1470)W*X»Y 
1450  PRINT 
1460  WRITE  < 15* 190) 

1470  FORMAT  4F15.7 

1480  IF  S=0  AND  Y=0  THEN  1530 

1490  GOTO  2140 

1500  STOP 


Figure  7*7  (continued) 
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1510  HEM 
1520  HEM 
1530  PRINT 

1540  HEM  THIS  PROGRAM  CALCULATES  THE  COMPLEX  ROOTS  OF  A FoiJPTH 

1550  PEN  ORDER  POLYNOMIAL  USING  THE  MODIFIED  BAIRSP'U  Ntr-iOr. 

1560  F;9=  1 

1570  39=1 

1530  A4=X 

1530  A3=W 

1600  A2»V 

1610  A1*U 

1620  A0*T 

1630  64*A4 

1640  B3*A3+R9*B4 

1650  B2=A£+R9*63+S9*B4 

1660  B1=AH>R9*B2+S9*B3 

1670  B0=A0+S9*B2 

1680  IF  ABS<B0>>lE-02  THEN  1700 

1690  IF  ABS<BlX1E-02  THEN  1730 

1700  S9=-A0/'B2 

1710  R9=-<A1+S9*B3>/B2 

1720  GOTO  1630 

1730  IF  <R9+2+4*S9><0  THEN  1730 
1 740  X 1 = < R9+SQR < P9t2+4*S9 ) ) /Z 
1 750  X2=< R9-SQR < F 9t2  +4*S9 > > /£ 

1760  X5=X6=0 

1770  GOTO  1810 

1780  Xl=X2=R9-2 

1790  X5*SQR<AB$<R9+2+4*S9)>/2 

1800  X6*-X5 

1310  IF  <B3t2-4*B4*B.2><0  THEN  I860 
1320  X3*<-B3*SGR  B3t2-4*34*B2> >/2 
1830  X4=<-B3-SGR' B3t2-4*34*B2> >/2 
1840  X7=X8=0 
1850  GOTO  1890 
I860  X3=X4*-B3/£ 

1870  X7*SQR< ABS<f3r2  .*64*62)  X2 
1880  XS*-X7 

1390  PRINT  "THE  EIGENVALUES  ARE" 

*900  PRINT 

1910  PRINT  " REAL  PART  IMAG  PART- 

1920  WRITE  <15*1990)X1»X5 

1930  WRITE  (15* 1990)X2*X6 

1940  WRITE  <15» 1990>X3»X7 

1950  WRITE  <15* 1990) X4» XS 

I960  PRINT 

1970  PRINT 

1980  WRITE  <15*190) 

1990  FORMAT  3X* F10. 5 * 5X» F10. 5 
2000  FIXED  4 
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2010  PRINT 

202O  PRINT  "MODE  I FREQUENCY  IS  "SGROil t2+X5t2> "RADIANS  PEP  SECOND" 
2030  PRINT  '*  PERIOD  IS  " 2#P I /SQR < X 1 t2+X5+2 ) ” SECONDS “ 

2040  PRINT  " DAMPING  RATIO  IS  “ -X 1 /SQR i X 1 t 2+X5*2 > 

2050  PRINT 

2060  PRINT  “MODE  I!  FREQUENCY  IS  “SQR<X3t2+X7t2) "RADIANS  PER  SECOND" 
2070  PRINT  " PERIOD  IS  " 2*P I /SQR < X3T2+X7 1 2 > " SECONDS " 

2080  PRINT  " DAMPING  RATIO  IS  B-X3/SQR<X3t2+X7t2) 

2090  PRINT 

2100  FORMAT  5F12.5 

2110  X9=X 

2120  Y9«Y 

2130  STOP 


Figure  7-7  (continued) 
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FEEDBACK  PARAMETER  IS  PITCH  RATE 
GA  I N = O 

FEEDBACK  IS  HOT  APPLIED 


******** 

THE  COEFFICIENTS 
f FROM  St<0> 

COEFFS  OF  G<  1 
1144. 7927540 

COEFFS  OF  G(  2 
-79. 9906392 

COEFFS  OF  G<  3 
0 . yOOOOOO 


*£4ifr,*,*jfrjr,fci£-ilr££iir-£ 

OF  THE  NUMERATOR  OF  THE  CLOSED-LOOP  TX  PH  ARE 
TO  St<3>> 

> 

309.71 38256  -0.111 6300  0 . OOOOOOO 

> 

-168.3869827  - 17360. 66 10600  -69.30000^0 


-0.3483321  -35.8048573  -26.0090000 


COEFFS  ‘I  j’  p > 

-0.  •*:  ‘ -35.8048573  -26.O090OOO  0.  OOOOOOO 


CQEFF 1 1. : L rS  .*»  THE  DENOMINATOR  ARE*. 

(FROM  1>  TO  ST<5>> 

0. ‘.OOT.joO  0.0779700  O.19329O0  19.5263700 

4 21 4 -yoy  1.0000000  0.000OOO0 

* * * ‘ * *************** 


THE  EIx  ■'  - *-9  ARE 


REA'  »A* 
-0.  0-J4S3 

-0.  01.1453 
-2. 10532 
-2. 10532 


IMAG  PART 
8.06308 
-O.06308 
3.87967 
-3.87967 


* 4 * * * * 


********* 


****** 


MODE  I FREQUENCY  IS  0.0632  RADIANS  PER  SECOND 
PERIOD  IS  99.3505  SECONDS 
DAMPING  RATIO  IS  0.0716 


MODE  II  FREQUENCY  IS 
PERIOD  IS 
DAMPING  RATIO  IS 


4.4141  RADIANS  PER  SECOND 

1.4234  SECONDS 

0.4770 


Figure  7-8 
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Figure  7-9  (a,'n,c,d)  shows  how  the  roots  vary 
as  the  feedback  sain  is  varied  and  as  different  para- 
meters are  fed  back.  Complex  conjugates  are  not  shown 
where  their  omission  would  be  confusing.  These  results 
are  tabulated  for  more  careful  scrutiny  in  Table  7-^» 

Figure  ?-10(a)  depicts  the  time  history  of  the 
''cRuer  F-89  aircraft  subject  to  an  initial  condition 
of  vertical  velocity  = 1 0 fus.  To  control  inputs  are 
postulated.  The  parameter  chosen  to  display  system 
response  is  pitch  angle;  it  was  chosen  because  of  its 
good  definition  of  short  period  and  phugoid  motion. 

Feedback  of  horizontal  velocity  error  can  sub- 
stantially increase  phugoid  damping’  before  the  short 
period  mode  is  much  affected.  However,  a further 
increase  in  the  gain  of  this  error  signal  will  drive 
the  short  period  to  instability.  The  sensitivity  of 
the  system  in  this  regime  is  shown  by  Figure  7-1 0(b); 
here,  with  ku  = 0.02,  the  short  period  motion  is  less 
than  critically  damped  and  oscillates  several  times 
before  dying  out.  The  damping  effect  on  the  phugoid 
is  readily  apparent. 

Feedback  of  vertical  velocity  (or  angle  of  attack), 
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f eedback 


(a)  ku 


Fimire  7-9 
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(d)  feedback 


Figure  7-9  (continued) 


133 


■3  ► 3 ^ O >s\  \/\  - 3 Pp  O ^ ^ H3  ^ ►-iJ  »~3  O %/\  ^-3  HI 


Feedback  parameter!  HORIZONTAL  VELOCITY  (u) 


0 

0.002 

0.02 

O 

• 

O 

T(sp) 

1 .42 

1 .44 

1 .60 

1.89 

"j  ( sp) 

0.477 

0.475 

0.444 

O.154 

T(ph) 

99.35 

17.86 

5.16 

2.76 

^ ( ph ) 

0.072 

0.097 

0.297 

0.702 

0.075 


.67 

-0.02 

2.55 

0.888 


1.52 

-0.10 

2.43 

0.980 


k parameters  VERTICAL  VELOCITY  (w) 


air. 

0 

m 

w 

.003 

( sp ) 

1.42 

1.77 

1.65 

1.32 

1.18 

0.73 

(sp) 

0.477 

1 

1 

\ 

■1 

1 

(ph) 

99.35 

165.17 

75.04 

68.11 

78.43 

89.15 

(ph) 

0.072 

1 

1 

0.108 

0.082 

0.072 

3 

- 

254.69 

112.48 

- 

- 

- 

3 

- 

-1 

- 

- 

- 

4 

- 

IO.54 

20.14 

9.56 

5.1° 

1.55 

4 

1 

1 

1 

-1 

• 1 

—1 

k par 

ameter  s 

! PITCH 

RATE  (q) 

air 

0 

0.01 

0.02 

0.05 

0.1 

0.25 

(sc) 

1 .42 

1 .41 

1 .40 

1.36 

1 . 51 

1.31 

(sp) 

0.477 

0.502 

0.526 

0.598 

0.710 

1 

(ph) 

99.35 

1 00.25 

1 01 .14 

103.75 

107.97 

1 1 9.69 

(ph) 

0.072 

0.072 

0.073 

0.074 

0.076 

0.084 

3 

- 

- 

- 

- 

- 

1.06 

3 

— 

* 

1 

k parameters  FITCH  ANGLE  (9) 
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if  the  rain  is  small,  will  increase  the  frequency  and 
damping  ratio  of  the  short  period  mode.  Since  the  angle 
of  attack  signal  is  l/66oth  of  the  vertical  velocity 
signal,  it  would  be  a better  parameter  from  which  to 
obtain  an  error  signal.  Increasing  the  gain  of  this 
feedback  will  cause  a third  mode  to  go  unstable;  however, 
its  period  is  so  great  that  no  problem  would  be  antici- 
pated in  controlling  it  manually.  A further  increase 
in  gain  will  cause  a fourth  mode  of  signif icantly  shorter 
period  to  go  unstable.  Since  its  period  decreases  ran  idly 
with  an  increase  in  gain,  it  may  be  assumed  that  values 
of  gain  in  this  region  are  to  be  avoided.  Figure  ?-'0 
(c)  shows  a time  history  of  pitch  angle  versus  time  with 
kw  = 0,0001  and  w0  = lOfps. 

Feedback  of  pitch  rate  has  almost  no  effect  upon 
the  phugoid  motion  but  appreciably  increases  the  short 
period  damping.  This  very  desirable  feature  makes  this 
parameter  useful  in  speed,  altitude,  and  attitude  hold 
modes  for  auto-pilots.  Figure  7-10(d)  shows  a time 
history  of  pitch  angle  against  time  with  kq  = 0.1  and 
wQ  = lOfps. 
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Firure 


Figure  7-10  (continued) 


Pitch  angle  feedback  is  also  an  excellent  method 
of  improving  longitudinal  response.  In  this  case,  the 
phugoid  damping  is  signif icantly  increased  at  the  expense 
of  short  period  damping.  If  short  period  damping  is  mar- 
ginal to  begin  with,  use  of  this  feedback  may  cause  an 
instability  in  that  mode.  Figure  7-1 0(e)  shows  the 
effects  of  kj  = 0.05  on  the  time  history  of  the  pitch 
angle  perturbations  with  w0  = lOfps. 

Use  of  both  pitch  rate  and  pitch  angle  feedback 
together  should  provide  an  ideal  system  with  maximum 
increase  in  damping  ratios  and  minimum  increase  in 
frequency.  Using  the  same  feedback  gains  as  above,  the 
response  is  shown  in  Figure  7-1 0(f). 

The  equations  of  state  for  each  of  the  above 
examples  were  obtained  by  entering  BASMAT  with  the 
modified  plant  matrix  ^A  - bkT}  and  multiplying  the 
resulting  state  transition  matrix  by  the  initial  condition 
vector  in  accordance  with  the  equation 

jc(t)  = $>(t)^c(0)  where  £(0)  = 

One  further  conclusion  was  reached  in  this  treatment 
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of  feedback!  the  two  parameters  which  are  mos^  affected 
in  phueoid  motion,  forward  velocity  and  pitch  angle,  were 
the  two  whose  feedback  most  affected  that  motion.  Their 
effect  on  the  short  period  was  minimal.  Likewise,  the 
two  parameters  which  are  relatively  constant  in  the 
phupoid  but  which  vary  considerably  in  the  short  period, 
anele  of  attack  and  pitch  rate,  had  little  effect  on 
the  phugoid  but  significantly  damped  the  short  period. 
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VIII 


The  Aerodynamic  D°s i gn  Data  Book  was  obtained  from 
TASA's  Lyndon  B.  Johnson  Space  Tenter  in  order  to  obtain, 
data  necessary  to  calculate  the  stability  derivatives  and 
to  acquire  the  block  diagram  of  the  longitudinal  control 
system  (automatic  mode)  of  the  srace  shuttle  orbiter. 

The  environment  chosen,  for  the  study  was  the  area  at 
which  the  Vach  number  was  " .05.  This  area  was  chosen 
because  the  aerodynamics  are  undergoing  raoid  change  here, 
a^d  it  would  provide  a testier  challenge  for  programming 
and  solution. 

~he  stability  derivatives  were  found  by  plotting 
the  force  and  moment  coefficients  with  respect  to  velocity 
(’7)  and  angle  of  attack  and  then  measuring  the  slope;  see 
Figures  P-1  through  8-6.  3y  choosing  points  close  together 
where  the  slope  was  rapidly  changing,  piecewise  linearity 
could  be  assumed  between  them.  The  results  are  shown  in 
Table  8-1.  The  derivatives  with  respect  to  61  were  also 
taken  from  the  publication;  the  graph  of  A CD<  -vs- 
is  shown,  as  an  example,  in  figure  8-7.  The  £t  derivatives 
are  listed  in  Table  8-2. 
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MriOO 


AOA  (dee-) 

Cl. 

5 

7.5 

0.0564 

I 0 

12.5 

1 5 

AOA  (dee-) 

£ D. 

5 

0.0046 

1 0 

0.0133 

15 

0.0234 

20 

AOA  (der) 

Cm. 

7 . 5 

-0.0030 

10 

-0.0043 

12.5 

-0.0036 

15 

-0.0043 

TJ  ( f p s ) 

Cl. 

920.1 

0.0 

1016.9 

7.23E-4 

IO65.3 

-1 .242-4 

1162.2 

-1  • 752-4 

U (fps) 

c.. 

968.5 

5.792-4 

1 016.9 

3 . 1 OH-4 

1065.3 

1 .24E-4 

1113.7 

-2.1  OS-5 

U (fps) 

92O.I 

-6 

1016.9 

-1.552-4 

1065.3 

0.0 

1162.2 

-1 .65E-4 

Table  3-1 
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0.0057 

0.0072 

0.0076 


0.0036 


20 


0.0*  03 
O.Qi 42 


St 

B9 

-25* 

-0.0013 

-20 

-0.0014 

-15 

-0.0012 

-10 

0.0024 

-5 

0.0050 

0 

0.0080 

— 

-30’ 

-0.0010 

-25 

-0.0016 

-20 

-0.0026 

-15 

-0.0042 

-10 

-0.0053 

0 

-0.0067 

20 

-O.OO89 

Table  3-2 
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In  cases  where  U was  the  independent  variable,  angle 
of  attack  was  assumed  to  be  "O0;  where  a.nele  of  attack 
was  the  independent  variable,  " was  1016.9  fps  (Kach  = 
1.05).  The  derivatives  were  cased  upon  an  anele  of 
attack  of  10°, 


The  point  in  the  desim  reentry  trajectory  where  the 
Mach  number  is  1.05  occurs  at  approximately  62000  reet  ir 
altitude.  Also,  at  this  point,  the  load  factor  n2  = 

1.09,  and  the  air  density  f - 2.0325H-4  slugs/ ft'. 


Other  constants  which  are  reauired  for  comoutations 


wins'  area  = 3 = 2690  ft^ 
wins'  K.A.C.  = c = 39*56  ft 
moment  of  inertia  = Iy  = 5.7330E+6  slue-ft^ 

gravity  = g = 32.073  ft/sec^  (corrected 
for  altitude  and  latitude) 


pitch  angle  = 0 = 9.0  dee 

mass  = m = 5339  slue's 


For  the  condition  of  Mach  = 1.05  and  anele  of  attack 
- <0°,  the  stability  derivatives  take  on  the  following 


values : 


- O.  0^4  > O 
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Entering  these  values  into  the  plant  matrix  (see  pare 


54)  yields 


-0.04410 

0.02849 

0.0 

-32.07 

-0.09845 

-0.01635 

1016.9 

0.0 

-4 , 683E-4 

-9.730E-6 

-0.09877 

0.0 

0.0 

0.0 

1 .0 

0.0 

f -0.44219  A 
-0.78489 
-0.01476 

, o » o 


Entering  rASl'AT  with  [aJ  yields  the  results  shown  ir. 
Figure  8-8.  Using  the  GRAPH  program  developed  in  Chapter 
VI  (see  page  71).  the  time  response  to  an  initial  condi- 
tion of  wQ  =10  fps  was  plotted.  The  results  are  shown 
in  Figure  8-9«  The  vehicle  is  clearly  unstable  in  the 
open-loop  conf iguration. 


A look  at  the  eigenvalues  printed  out  in  Figure  8-8 
showed  that  the  phugoid  had  disappeared,  being  replaced 
by  two  real  roots,  one  of  which  was  positive.  Although 
the  negative  (stable)  root  was  larger,  the  instability 
became  dominant  very  quickly. 

The  eigenvalues  were  not  widely  separated,  having 
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Figure  8-8 
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Pigure  8-8  (continued) 


continued ) 


periods  of  75. 1 , 61.2,  and  39*7  seconds,  respectively, 
so  the  short  period  was  effectively  missing  as  well. 
Although  its  theoretical  damping  ratio  was  0.^337,  it 
was  completely  supressed  by  the  instability  of  the 
positive  real  root. 

It  was  obvious  that  some  form  of  feedback  network 
would  be  necessary  so  that  the  stability  of  the  vehicle 
would  be  more  in  line  with  that  of  conventional  aircraft. 
And,  since  the  orbiter  was  designed  to  fly  at  the  very 
extremes  of  atmospheric  conditions,  the  system  would 
necessarily  depend  upon  many  parameters,  among  them 
dynamic  pressure,  Mach  number,  true  airspeed,  load  factor 
(which  is  a function  of  Ct ) , pitch  angle  and  rate,  angle 
of  attack,  and  elevor.  deflection  angle.  Other  factors 
could  have  been  used,  e.g.,  vertical  velocity,  descent 
angle,  to  name  a few,  but  the  ones  chosen  appeared  to  be 
those  most  easily  obtainable  and  from  which  other 
necessary  information  could  be  computed.  The  end  result 
of  the  engineering  design  is  shown  in  Figure  8-10. 

In  order  to  solve  the  equations  of  motion  as  they 
were  modified  by  the  various  gains,  feedbacks,  lead-lag 
elements,  and  limits,  the  Continuous  System  Modeling 
Program  (version  III),  as  incorporated  on  the  IBM360, 
was  utilized. 
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Firure  8-10 


The  system  variables  were  assigned  names  as  follows: 

Nz  command  = C-CMD 
true  airspeed  = U1 
Nz  ~ g = NZ 

$ = THETAl 

horizontal  velocity 

perturbation  = XI 

vertical  velocity 

perturbation  = X2 

pitch  rate  perturbation  = X3 
pitch  angle  perturbation  = X4 

elevon  deflection 

perturbation  = X8 

The  input  signal  to  the  first  lead-lag  element  was 
called  X5I':  and  its  output  called  X5 • The  input  to  the 
second  lead-lag  element  was  called  X6lN  and  its  outpu*  Xc. 

The  signal  entering  the  St  LIMIT  was  called  XQFR~’. 
and  the  limited  sienal  X3FRE2:  the  sipr.al  enterin'  '.te 
<Se  LIMIT  was  called  X3FRE3  and  its  output  called  XQDC7 : 
the  integration  of  X8DCT  was  X3.  The  signal  XR  was  fed 
back  through  a lap-  element  whose  output  was  labeled  X?. 

The  CSThP  III  program  is  divided  into  three  segments, 
INITIAL,  DYNAMIC,  and  TERMINAL.  The  INITIAL  part  is  for 
defining  parameters  whose  value  will  not  change  for  the 
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duration  of  the  program.  The  DYNAMIC  segment  contains 
the  equations  of  motion  and  other  computations.  The 
TERMINAL  section  (optional)  contains  any  final  operations 
which  are  desired.  For  this  analysis,  the  following 
terms  were  assumed  constants 

gravity  = GRVTY  = 32.073  ft/sec2 
Xq  = ZQ  = 0.0 

density  = RHO  = 2.0325E-L  slue/ft3 
Cv  = CMADOT  = -2.0 

The  initial  values  of  two  parameters  were  listed: 
pitch  angle  = THETA  = 0.1571  rad 
true  airspeed  = UO  = 101 6. 9 fps 
The  other  constants  listed  in  the  program  (see  Figure 
9-11)  are  dimensional  parameters  or  mass  properties  of 
the  orbiter.  AREA  refers  to  wetted  wing  area;  I3UBY 
(Iv)  is  the  mass  moment  of  inertia  about  the  pitch  axis; 
CHORD  is  the  reference  wine  chord  leneth;  MASS  is  the 
mass  of  the  vehicle  with  the  payload  out. 

The  FUNCTION  statements  enable  the  proeram  to  find 
a dependent  variable  knowing  the  value  of  the  independent 
variable,  much  like  reading  a graph.  The  interpolation 
between  the  points  listed  in  the  statement  may  be  linear 
or  quadratic.  Should  an  independent  variable  value  fall 
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outside  the  ran?e  of  points  listed,  the  program  will 
extrapolate  a value  based  on  the  slope  of  the  function 
at  the  nearest  point  listed. 

The  FUNCTION  statements  are  entered  with  an  AFC-EN 
command  (for  linear  interpolation).  Since  linear  inter- 
polation was  used  in  this  analysis,  it  was  necessary  to 
insure  that  the  dependent  variable  varied  linearly 
between  the  points  chosen.  The  functions  were 
FI  i CMDE  - Cm* 


F2  : 

CODE  = 

Cp* 

F3: 

CEDE  = 

F4: 

GCKD  = 

F5  * 

CLA  = 

Cw. 

F6  i 

CD  = 

Go 

F7 « 

CLU  = 

FSi 

CL  = 

c. 

F9  i 

CDA  = 

Co. 

F1  Oi 

CDU  = 

FI 1 : 

f'MTT  — 
vl’iU  — 

FI  2: 

= 

Cm 

FT  3 i 

CMA  = 

Cm. 

FI  4 i 

CI'/’Q  = 

C*«y 
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These  coefficients  were  in  turn  used  to  compute  the 
stability  derivatives.  As  the  values  of  the  independent 
variables  changed,  so  did  the  dependent  variables  and  so 
also  did  the  stability  derivatives.  In  this  manner,  the 
stability  derivatives  were  modified  due  to  changing 
flight  conditions,  and  thus  an  important  facet  of  the 
model  was  maintained.  The  stability  derivatives  were 
then  used  to  write  the  equations  of  motion  (the  XI DOT 
through  X^DOT  equations). 

The  modeling  of  the  system  was  completed  by  writing 
four  other  equations  as  the  block  diagram  signified. 

X5  and  X6  were  written  simply  as  lead-lag  (LSDLAS) 
functions  of  X5 I"  and  X6IX,  respectively,  where  the 
parameters  listed  were  the  coefficients  of  s in  the 
numerator  and  denominator. 


The  term  X7  was  found  by  writing  the  lag 
from  which  it  came  as  a differential  equation 
output  divided  by  input  equals  the  transfer  f 


X? 

X^ 


-Lil 


s + 1.5 


(X7)s  + l.5(X7)  = l.5(X8) 


compensation 
. Since 
unction, 


Taking  the  inverse  Laplace  transform, 
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X*7  = 1.5(X9)  - l.5(X7) 
and  X7  = J X*7 

Writing  this  in  the  CS"T  III  language, 

X7  = ITT  GRL ( 0 . 0 , X 7D0T ) 

where  the  0.0  refered  to  the  initial  condition. 

The  equation  to  find  X3  was  somewhat  more  complex. 

The  signal  was  first  the  sum  of  GD3*X6  + X7  = X3F3E1  ? 
this  was  then  limited  to  -35/+200,  compared  to  XS,  and 
the  difference  multiplied  by  20  (here,  in  radians). 

After  limiting  this  signal  to  -20/+200  per  seco"d,  it  was 
integrated  to  find  X8. 

The  eight  equations  were  solved  simultaneously  using 
the  STIFF  method  of  integration,  and,  after  converting 
some  of  the  parameters  into  degrees,  the  results  were 
printed . 

It  may  be  seen  from  the  block  diagram  that  the  forcing 
function  is  the  command  load  factor,  GCKD.  This  value, 
plotted  as  FUNCTION  F^ , was  available  from  the  Aero- 
dynamic D^s i m Data  ~ ook  and  was  programmed  accordingly. 
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The  value  of  CMADOT  was  not  available  from  the 
above  publication.  To  arrive  at  the  value  of  -2.0, 
several  aircraft  were  studied,  considering  wing  plan- 
form,  aircraft  size,  airspeed,  and  altitude.  The  larger 
aircraft  (C-47  and  DC-8)  had  values  of  -10  and  -6.8, 


respectively,  the  DC-8 

figure 

being  a 

t Mach 

0.83  and 

33000  feet.  The  A-4D, 

having 

a delta 

plan.-1’ 

orm,  had  a 

Cyk  of  -2.1  at  15000  f 

e e t and 

-1.4  at 

35000 

f QO  + « 

Considering  the  size  of  the  orbiter  and  its  altitude 
for  this  study,  the  value  of  -2.0  was  deemed  appropriate. 

Several  equations  to  compute  the  load  factor  on  the 
vehicle  were  considered.  The  one  chosen  was  derived  from 
the  lift  equations 

nz  = ^ where  L = Ci,q3 

so,  nz  = C^aS/mg 

= CL*QHAT*AREA/(MASS*GRVTY) 

This  equation  was  chosen  as  it  was  most  independent  of 
the  equations  of  motion  and  would  therefore  amplify  to 
a less  extent  any  errors  due  to  integration  round-off  or 
to  other  sources  during  the  solution  of  the  simultaneous 
equations . 


169 


The  results  are  shown  in  Figure  8-12  (short  period) 
and  Firure  8-13  (phuvoid).  Inasmuch  as  both  modes  were 
evidently  stable,  the  steady  state  magnitudes  were 
extremely  high.  Several  causes  must  be  considered. 

One,  the  equations  of  motion  were  written  as  per- 
turbed motions.  This  was  necessary  in  order  to  justify 
their  linearization.  The  parameter  values  were  treated 
as  perturbations  throughout  the  program,  being  added  to 
the  larre  scale  motion  where  necessary  for  computation 
(as  in  the  case  of  THETA1  and  U1 ) and  for  output  plots. 
Their  initial  conditions  were  collectively  zero.  The 
initial  conditions  for  the  other  parameters,  X5  through 
X?,  were  not  zero,  however,  and  a method  for  incorpora- 
ting (or  even  deriving)  their  initial  value  was  not 
available.  The  initial  condition  for  X8  was  available 
from  the  Data  Book  and  was  so  used.  It  is  not  known 
if  the  lack  of  initial  conditions  for  these  parameters 
caused  the  large  excursions  in  the  response?  it  is 
doubtful  that  this  was  the  case  because  the  steady  state 
values  should  still  have  been  reasonable;  they  were  not. 

Two,  the  simulation  was  be.sun  at  some  point  other 
than  t=0  in  the  trajectory.  That  is,  the  vehicle  had 
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Fipure  H-1 2 (continued) 
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Figure  8-12  (continued) 
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Fi^rure  8-12  (continued) 
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Figure  8-12  (continued) 
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Figure  8-12  (continued) 
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Figure  8-12  ( co-i  Linued ) 
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Figure  8-12  (continued) 
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Figure  8-1  3 
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Figure  8-13  (continued) 
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Figure  8-13  (continued) 
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Figure  8-13  (continued) 
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Figure  8-13  (continued) 
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Figure  8-13  (continued) 
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Fifnire  8—1 3 (continued) 


0*0 


Figure  8-13  (continued) 


already  entered  the  atmosphere,  was  flying  aerodynamically, 
being  commanded,  responding,  etc.,  in  short,  was  in  a 
dynamic  environment  when  this  program  started.  Unless 
the  parameters  existing  at  t0  for  this  simulation  can  be 
precisely  ascertained,  it  is  doubtful  that  meaningful 
results  can  be  expected. 

To  this  end,  angle  of  attack  corresponding  to  the 
commanded  load  factor  at  t0'  was  calculated t 

nz  * 1.09  » Ci,(ifV2)s/mg 

Solving  for  yielded  0.722.  Since  the  angle  of  attack 
at  1.0  "g"  was  10°  with  a Cl  of* 0.544,  the  angle  of 
attack  at  1.09  "g"  was  found  to  be  10°( 0.722/0. 544) 
or  13.27°  (0.2316  rad). 

An  initial  pitch  rate  would  affect  the  X2  coefficient 
in  the  X1D0T  equation  and  the  XI  coefficient  in  the  X2D0T 
equation,  so  it  was  computed  using  nz  a Q0U0/g  +1.0 
and  found  to  be  0.0344  rad/sec. 

Both  of  these  modifications  were  made  to  the  program 
with  negligible  effect  on  the  results. 

To  analyze  the  results  as  they  standi 

The  angle  of  attack  first  exceeded  the  limits  of  the 
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FUNCTION  statements  (requiring  interpolation  of  the  force 
coefficients  well  outside  the  linear  range  of  the  functions) 
at  approximately  2 seconds*  oscillated  in  and  out  of  the 
range,  and  left  the  range  for  good  at  approximately  10 
seconds*  All  stability  derivatives  dependent  upon  angle 
of  attack  must  be  considered  suspect  after  this  time* 

The  angle  of  attack  was  computed  by  dividing  the  vertical 
velocity  perturbation  by  the  horizontal  velocity  and  adding 
this  to  the  initial  angle  of  attack. 

The  vertical  velocity  response  (positive  is  downward) 
reached  an  extremely  high  value  very  quickly,  and  if  applied 
to  altitude,  would  have  the  orbiter  impacting  the  ground 
before  the  simulation  period  was  half  over.  Such  a rate 
of  descent  would  also  cause  a rapid  change  in  air  density, 
which  was  originally  assumed  constant,  leading  to  very 
different  values  of  those  stability  derivatives  proportional 
to  RHO.  Tv  assumption  of  a constant  RHO  should  be  valid 
”nr  the  •*  .'’ee  minute  duration  of  the  simulation. 

me  pitch  rate  magnitudes  are  not  unusual,  perhaps  a 
bit  high,  but  the  steady  state  pitch  angle  is  much  too 
high.  The  steady  state  values  of  58°  pitch  angle  and  71° 
angle  of  attack  indicate  a glide  angle  of  13°,  also  much 
too  high. 
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The  elevon  deflection,  going  to  the  trailing  edge  up 
limit  in  just  over  ten  seconds,  indicated  that  the  vehicle 
sensed  an  insufficient  nose  up  pitching  moment  with  a 
GCMD  of  1.09  and  was  trying  vainly  to  create  that  moment. 
The  result  of  the  full  trailing  up  deflection  was  instead 
a stalling  angle  of  attack  and  a high  rate  of  descent. 

The  actual  load  factor  experienced  by  the  vehicle  in 
the  simulation  was  only  slightly  higher  than  expected*  the 
laree  elevon  deflection  was  undoubtedly  the  cause  of  the 
initial  rapid  rise  in  this  parameter.  However,  since  the 
program  computed  nz  as  a function  of  p , the  steady  state 
value  was  meaningless. 

Taking  the  above  observations  into  consideration,  the 
program  appeared  to  be  consistent.  Many  possible  reasons 
for  the  magnitude  discrepancies  were  examined,  but  no 
changes  in  the  program  made  more  than  a modicum  of  change 
in  the  results.  It  appeared  that  the  error  was  basic  in 
nature  and  not  related  to  the  model  as  presented  or  to 
the  CSMP  proerram.  The  most  likely  cause  is,  as  was 
mentioned  earlier,  the  fact  that  the  simulation  is  not 
commencing  at  t=0. 
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